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Recently,  variational  and  PDE-based  methods  have  become  very  popular  in 
the  processing  and  analysis  of  medical  images.  In  this  presentation,  we  apply  these 
methods  to  three  key  topics  in  medical  image  analysis  namely  i)the  correspondence 
problem,  ii)fuzzy  segmentation  And  iii)  registration. 

The  Correspondence  Problem  for  curves  has  many  applications  to  medical 
imaging,  including  motion-tracking,  morphometries  and  neural  growth.  The  general 
problem  is  to  find  a continuous,  non-rigid  correspondence  between  two  curves  such 
that  a shape-dissimilarity  measure  is  minimized.  Current  methods  solve  this  problem 
to  its  entirety  only  for  parameterized  curves.  In  this  work,  we  look  at  the  correspon- 
dence problem  for  implicitly  represented  curves.  Given  two  signed  distance  functions, 
we  search  for  a diffeomorphism  between  their  zero-level  sets  that  minimizes  a shape 
dissmilarity  measure.  The  diffeomorphisms  are  generated  as  flows  of  certain  tangen- 
tial vector  fields,  and  curve-normals  are  chosen  as  the  similarity  criterion.  We  also 
show  how  this  model  can  be  extended  to  compare  curves  of  different  topologies.  We 
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have  tested  our  model  on  synthetic  and  ultrasound  cardiac  data,  and  have  obtained 
good  results. 

Segmentation  of  images  is  a widely  researched  field  and  is  unarguably  the 
first  step  to  any  problem  in  medical  image  analysis.  In  this  work,  we  use  fuzzy 
classification  and  active  contours  in  a single  variational  framework,  allowing  the  use 
of  tools  from  both  deformable  geometry  and  clustering,  and  hence  giving  a useful 
technique  for  unsupervised  segmentation.  The  model  was  tested  on  synthetic  and 
MRI  brain  data,  with  promising  results. 

For  techniques  like  functional  MR  imaging  that  work  with  sequences  of  im- 
ages, proper  registration  or  spatial  alignment  of  the  images  is  essential.  A usual 
registration  technique  is  to  segment  a key  feature  in  each  of  the  images  and  then 
align  the  images  with  respect  to  this  feature.  But  due  to  low  contrast  and  poor 
resolution  of  fMR  images,  usual  techniques  that  use  only  image  information  for  the 
segmentation  step,  give  inaccurate  registration  results.  In  our  work,  we  incorporate 
prior  shape  information  into  the  segmentation  step,  within  a variational  framework 
and  have  obtained  good  registration  results. 
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CHAPTER  1 
INTRODUCTION 


Rapid  advances  in  Medical  imaging  technologies  have  significantly  improved 
the  quality  of  a patient’s  medical  care  today.  While  planar  X-ray  imaging  was  the 
only  radiological  imaging  method  in  the  early  part  of  the  last  century,  several  mod- 
ern imaging  modalities  are  in  practice  today  to  acquire  anatomical,  physiological, 
metabolic  and  functional  information  from  the  human  body.  This  has  allowed  physi- 
cians to  make  increasingly  accurate  diagnosis  and  render  precise  and  measured  modes 
of  treatment.  Current  applications  of  imaging  technologies  include  laboratory  medi- 
cine, surgery,  radiation  therapy,  nuclear  medicine,  and  diagnostic  radiology. 

The  commonly  used  medical  imaging  modalities  capable  of  producing  multidi- 
mensional images  are:  X-ray  computed  tomography  (X-ray  CT),  magnetic  resonance 
imaging  (MRI),  single  photon  emission  computed  tomography  (SPECT),  positron 
emission  tomography  (PET)  and  ultrasound.  The  above  methods  make  use  of  sophis- 
ticated instrumentation  and  equipment  using  high-speed  electronics  and  computers 
for  acquisition,  reconstruction  and  display  of  digital  medical  images.  These  digital 
images  of  physiological  structures  can  then  be  processed  and  manipulated  to  reveal 
hidden  characteristic  diagnostic  features.  Further,  these  features  of  interest  can  be 
quantified  and  analyzed  to  understand  their  behavior  to  help  with  a diagnosis  or  to 
evaluate  a treatment  protocol. 

In  general,  all  medical  imaging  modalities  employ  an  external,  internal  or  a 
combination  of  energy  sources  for  image  aquisition.  Imaging  methods  like  X-ray  CT, 
usually  make  use  of  ionized  radiation  such  as  X-rays  as  an  external  energy  source  and 
are  primarily  used  for  anatomical  imaging.  Such  anatomical  imaging  modalities  are 
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based  on  the  attenuation  coefficient  of  radiation  passing  through  the  body.  Another 
example  of  external  energy  source  based  imaging  is  ultrasound  or  acoustic  imaging. 

Nuclear  Medicine  imaging  modalities  such  as  SPECT  and  PET  use  an  inter- 
nal energy  source  through  an  emission  process  to  image  the  human  body.  These 
modalities  provide  useful  metabolic  information  about  the  physiological  functions  of 
the  organs.  For  emission  imaging,  radioactive  pharmaceuticals  are  injected  into  the 
body  to  interact  with  selected  body  matter  or  tissue  to  form  an  internal  source  of 
radioactive  energy  that  is  used  for  imaging. 

Further,  a clever  combination  of  external  stimulation  on  internal  energy  sources 
can  be  used  in  medical  imaging  to  acquire  more  accurate  information  about  tissue 
material  and  physiological  responses  and  functions.  MRI  uses  external  magnetic  en- 
ergy to  stimulate  selected  atomic  nuclei  such  as  hydrogen  protons.  The  excited  nuclei 
become  the  internal  source  of  energy  to  provide  electromagnetic  signals  for  imaging 
through  the  process  of  relaxation.  The  relaxation  process  is  controlled  by  the  bio- 
logical parameters  Tj  and  T2.  These  parameters  are  tissue  dependent,  introducing 
the  possibility  to  separate  different  tissue  types  in  the  human  body.  MRI  of  the 
human  body  provides  high-resolution  images  of  the  human  body  with  excellent  soft- 
tissue  characterization  capabilities.  Recent  advances  in  MRI  have  led  to  perfusion 
and  functional  imaging  aspects  of  human  tissue  and  organs  [77] . 

Image  processing  and  analysis  methods  are  aimed  at  the  enhancement  of  di- 
agnostic information  to  improve  interpretation  of  medical  images.  Nowadays,  the 
size  and  number  of  these  medical  images  have  necessitated  the  use  of  computers 
to  automate  their  quantitative  processing  and  analysis.  Methods  of  processing  and 
analysis  include  tasks  such  as  reconstruction,  noise-removal  and  image-enhancement, 
segmentation,  shape  representation  and  matching,  and  motion-tracking. 
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Image  analysis  methods  constitute  an  eclectic  collection  of  techniques  derived 
from  many  different  theoretical  standpoints  such  as  signal  processing[53],  mathe- 
matical morphology [76],  stochastic  approaches ([27,  48])  and  recently  variational  and 
partial  differential  equation  methods([l,  74]). 

In  this  presentation,  we  apply  variational  and  partial  differential  equation(PDE) 
methods  to  three  key  topics  in  medical  image  analysis,  namely  shape  matching,  seg- 
mentation and  registration. 

Variational  and  PDE  methods  have  been  widely  used  in  image  processing  in 
the  past  few  years.  Examples  include  continuous  mathematical  morphology,  shape 
from  shading,  invariant  shape  analysis,  image  segmentation,  tracking,  image  restora- 
tion, stereo,  contrast  enhancement  and  inpainting.  Using  variational  and  PDE  meth- 
ods in  image  analysis  leads  to  modeling  images  in  a continuous  domain.  This  simpli- 
fies the  formalism,  which  becomes  grid  independent  and  isotropic.  Conversely,  when 
the  image  is  represented  as  a continuous  signal,  PDEs  can  be  seen  as  the  iteration  of 
local  filters  with  an  infinitesimal  neighborhood.  This  interpretation  of  PDE’s  allows 
one  to  unify  and  classify  a number  of  the  known  iterated  filters,  as  well  as  to  derive 
new  ones.  Another  important  advantage  of  the  variational  and  PDE  approach  is  the 
possibility  of  achieving  high  speed,  accuracy,  and  stability  with  the  help  of  the  exten- 
sive results  of  numerical  PDE  methods.  Further,  this  area  has  a quite  unique  level  of 
formal  analysis,  giving  the  possibility  to  provide  not  only  successful  algorithms  but 
also  useful  theoretical  results  like  existence  and  uniqueness  of  solutions. 

In  the  following  chapters,  we  discuss  the  motivation,  formulation  and  numer- 
ical experiments  for  the  following  topics: 

1.  Shape-Matching  for  implicit  curves:  Finding  non-rigid  shape-based  correspon- 
dences between  implicitly  represented  curves. 
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2.  Fuzzy  Segmentation : A PDE  based  framework  that  combines  fuzzy  classifica- 
tion and  active  contours  to  give  a simple,  unsupervised  segmentation  technique, 
with  geometric  modelling  capabilities. 

3.  Functional  MRI  Registration : Joint  segmentation  and  registration  of  2D  fMR 
brain  images  using  the  corpus  callosum  as  a prior  shape. 

1.1  The  Non-Rigid  Correspondence  Problem  for  Implicitly  Represented  Curves 

Finding  non-rigid,  shape-based  correspondences  between  curves  has  many  ap- 
plications to  medical  imaging,  including  shape-matching,  motion-tracking,  morpho- 
metries (study  of  shape  variations)  and  neural  growth.  The  general  problem  is  to 
find  a meaningful  non-rigid  correspondence  between  two  curves  such  that  a shape- 
dissimilarity  measure  is  minimized.  The  problem  by  definition  assumes  the  given 
curves  to  be  noiseless  and  attempts  to  establish  a point-correspondence  between  the 
curves.  Current  methods  solve  the  problem  in  its  entirety  only  for  parametrized 
curves.  These  ideas  have  also  been  extended  to  find  shape-based  correspondences 
between  parametrized  surfaces.  However,  it  is  nontrivial  to  extend  these  methods 
to  compare  curves  of  arbitrary  topology.  In  the  medical  setting,  one  often  deals 
with  scenarios  that  require  comparing/tracking  curves  of  arbitrary  topology  (for  e.g. 
matching  two  spinal  columns  or  tracking  a cell  before  and  after  division).  This  mo- 
tivates us  to  consider  the  general  problem  of  finding  shape-based  correspondences 
between  curves  of  arbitrary  topology.  We  assume  that  the  curves  in  question  are 
implicitly  represented  as  0-level  sets  of  certain  2D  functions. 

A lot  of  work  has  gone  into  developing  techniques  that  handle  implicit  represen- 
tations in  vision  problems.  Recently  many  research  groups  have  shown  interest  in 
solving  the  curve-matching  problem  for  implicitly  represented  curves.  In  general, 
such  methods  use  a search  space  of  spatial  transformations(to  deform  the  embedded 
curves)  that  minimizes  the  Euclidean  distance  between  the  curves.  This  approach 
is  suitable  for  problems  such  as  image  registration  where  feature-dissimilarity  in  the 
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images  is  used  to  drive  the  required  transformation  between  the  images.  An  actual 
point-correspondence  between  the  features  is  not  required.  However  this  technique 
seems  to  be  an  overkill  for  curve-matching,  since  one  searches  a larger  space  of  spatial 
transformations  while  only  transformations  between  the  curves  need  to  be  considered. 
Further  these  methods  do  not  guarantee  point-point  correspondence  which  is  essen- 
tial in  many  applications  (e.g., motion  tracking).  One  of  the  novelties  of  our  work  is 
that  we  solve  the  correspondence  problem  for  implicitly  represented  curves  by  directly 
computing  transformations  between  the  curves.  Thus,  unlike  current  methods  which 
are  suitable  only  for  curve-matching,  our  work  is  applicable  both  to  curve-matching 
and  motion-tracking  of  implicitly  represented  curves. 

Given  the  functions  <f>i  : (fix  C !R2)  ->■  and  <f>2  : (^2  C 3fJ2)  SR  whose  0- 
level  curves  Ci  and  C2  we  want  to  match,  we  search  for  a diffeomorphism  that  aligns 
the  normal  vector  fields  of  C;.  The  problem  is  posed  in  an  variational  framework. 
We  model  our  search  space  of  correspondences  as  bimorphisms(originally  defined  in 
Tagare’s  work  [81]),  which  are  curves  on  the  torus(Ci  xC2), satisfying  certain  technical 
conditions.  Our  objective  functional  is  defined  as  an  arclength  integral  over  the  space 
of  bimorphisms  and  is  symmetric  with  respect  to  the  curves.  The  minimum  energy 
taken  by  this  objective  functional  defines  a notion  of  distance  between  the  curves 
C;  and  is  zero  when  the  curves  differ  only  by  scaling,  local  tangential  stretching 
and  rigid  motions.  A level-set  representation  for  the  functional  is  used  to  provide  a 
convenient  Cartesian  domain  in  which  to  perform  computations.  The  search  space  of 
bimorphisms  is  characterized  as  flows  of  certain  tangential  vector-fields.  This  setting 
allows  modelling  of  large  diffeomorphic  deformations  directly  between  the  curves. 

We  also  extend  the  above  framework  to  track  curves  of  different  topology  and 
see  a possibility  of  extending  it  to  3 dimensions  (work  in  progress).  The  model  was 
tested  on  synthetic  and  MRI  cardiac  data,  with  good  results. 
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1.2  A PDE-Based  Method  for  Fuzzy  Segmentation 

The  problem  of  dividing  the  image  domain  into  parts  which  are  homogeneous 
with  respect  to  some  image-feature(s)  is  termed  image-segmentation.  It  is  unarguably 
the  first  step  and  is  a very  important  aspect  of  image-analysis.  Segmentation  can  be 
broadly  classified  into  boundary-based  and  region-based  methods.  Boundary-based 
methods  rely  on  finding  pixels  with  high  image-gradients  and  then  grouping  them 
into  object  boundaries([10,  43,  52]).  Region-based  methods  classify  pixels  into  regions 
that  are  homogeneous  under  image-characteristics  ([64,  84,  97]). 

For  most  of  the  medical  imaging  modalities,  segmentation  is  a hard  problem 
to  solve  due  to  the  presence  of  noise,  artifacts,  low-contrast  and/or  -resolution.  There 
are  many  segmentation  approaches  for  medical  images  proposed  in  the  literature  that 
vary  depending  on  the  application,  noise-level  and  imaging  modality.  Of  these,  active 
contours  and  clustering-based  methods  are  commonly  used. 

Clustering  methods  [5]  are  pattern  recognition  techniques  that  classify  im- 
age pixels  into  a fixed  set  of  classes  based  on  the  feature-information  present  at  the 
pixel(e.g.  image  intensity).  These  methods  are  unsupervised  classifier  methods  in 
the  sense  that  no  prior  training  data  are  necessary  for  segmentation.  Clustering 
algorithms  such  as  K-means  are  popular  due  to  their  ease  of  implementation  and 
computational  speed.  But  the  lack  of  spatial  or  shape  constraints  on  the  segmented 
regions  makes  these  methods  susceptible  to  noise  and  intensity  inhomogeneities.  Ac- 
tive contours([93,  87])  are  physically  motivated,  model-based  techniques  that  can  be 
used  to  capture  closed,  Lipchitz  boundaries  of  regions  by  deforming  a parametric 
curve  under  the  influence  of  internal  and  external  forces.  The  internal  forces  con- 
strain the  deforming  contour  to  be  smooth  and  the  external  forces  pull  the  contour 
towards  region  boundaries. 

Active  contours  have  the  advantage  that  they  offer  a coherent  and  consistent  math- 
ematical description  and  are  robust  to  noise,  spurious  edges  and  boundary  gaps  due 
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to  their  incorporation  of  a smoothness  constraint.  They  can  be  extended  to  have 
topological  freedom  by  using  implicit  representations  [65]  instead  of  explicit  para* 
metrization.  But  active  contour  models  have  the  disadvantage  that  the  solution  is 
dependent  on  initial  contour  placement  and  convergence  can  be  slow  especially  if  the 
initial  contour  is  placed  far  away  from  the  image  boundaries. 

State-of-the-art  segmentation  techniques  combine  probability  based  or  clustering 
based  techniques  with  active  contours[80].  This  allows  use  of  region-based  statis- 
tics in  the  external  force  term  of  the  active  contour,  resulting  in  a fast  and  stable 
convergence  to  a global  minimizer. 

In  medical  imaging  modalities  such  as  MRI  and  CT,  the  sampling  process  results 
in  partial-volume  effects  leading  to  structural  ambiguities.  Partial-volume  effects  are 
artifacts  that  occur  when  multiple  tissue  types  contribute  to  a single  pixel,  resulting 
in  a blurring  of  intensity  across  boundaries.  A common  way  to  handle  this  problem  is 
to  allow  for  the  segmented  regions  or  classes  to  overlap,  termed  as  soft  segmentation. 
Fuzzy  classification  method  (FCM)  is  a well  known  fuzzy  counterpart  to  K-means. 

In  this  chapter,  we  use  a variational  approach  to  automatically  segment  med- 
ical images  into  a fixed  number  of  classes.  These  classes  form  a partition  of  the 
image-domain  into  connected  regions.  They  are  defined  such  that  intensity  values 
within  a class  are  close  to  a constant. 

We  have  used  FCM  and  active  contours  in  a single  variational  framework. 
This  allows  the  use  of  tools  from  both  deformable  geometry  and  fuzzy  clustering, 
thus  giving  us  a useful  technique  for  unsupervised  segmentation.  The  model  has 
shown  a high  robustness  to  noise  and  is  able  to  recover  complex  structures.  We  have 
tested  our  model  on  synthetic  and  T1  MRI  brain  data  and  compare  the  model’s 
performance  with  some  of  the  current  segmentation  methods. 
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1.3  Shape-Based  Segmentation  and  Registration  of  fMR  Images 

For  techniques  like  fMR  imaging  that  work  with  sequences  of  images,  proper 
registration  or  spatial  alignment  of  the  images  is  essential.  A usual  registration 
technique  is  to  segment  a key  feature  in  each  of  the  images  and  then  align  the  images 
with  respect  to  this  feature.  But  due  to  real-time  acquisition,  the  contrast  and 
resolution  of  the  images  are  very  low.  So,  usual  segmentation  techniques  that  use 
only  image  information  result  in  inaccurate  registration. 

Fortunately,  in  many  medical  imaging  problems,  prior  information  (e.g.  shape) 
of  the  feature  to  be  detected  is  available.  In  our  work,  we  have  combined  this  infor- 
mation with  an  edge-detection  term  in  a variational  framework  and  have  obtained 
good  segmentation  results.  In  addition  to  this,  the  rigid  transformation  that  aligns 
the  prior-shape  with  the  segmented  feature  is  estimated,  thus  achieving  simultaneous 
segmentation  and  registration.  We  have  tested  the  model  on  sequences  of  2D  fMR 
brain  and  ultrasound  cardiac  images,  obtaining  good  results. 


CHAPTER  2 

THE  CORRESPONDENCE  PROBLEM  FOR  IMPLICITLY-DEFINED  CURVES 

We  present  a novel  variational  model  to  find  shape-based  correspondences 
between  two  sets  of  level  curves.  While  the  usual  correspondence  techniques  work 
with  parametrized  curves,  we  use  a level-set  formulation  that  enables  us  to  handle 
curves  with  arbitrary  topology.  Further,  the  framework  can  be  extended  to  three 
dimensions.  Given  the  functions  dq  : (f^  C 5R2)  -4  3?  and  <f>2  : (D2  Q K2)  — > 5?  whose 
0-level  curves  we  want  to  match,  we  search  for  a diffeomorphism  that  minimizes  the 
rate  of  change  of  the  difference  in  tangential  orientation  of  the  zero-level  sets.  To 
make  the  formulation  symmetric  and  to  simplify  computations,  we  map  the  domains 
of  the  level-set  functions  <Iq  to  a common  domain  Q by  initial  diffeomorphisms  that 
are  chosen  arbitrarily.  We  then  search  for  diffeomorphisms  from  Q to  itself,  generating 
them  by  flows  of  certain  vector  fields  on  Q.  The  resulting  correspondences  are  scale- 
and  rotation-invariant  with  respect  to  the  curves.  We  show  how  this  model  can  be 
used  as  a basis  to  compare  curves  of  different  topology.  We  have  tested  model  was 
tested  on  synthetic  and  MRI  cardiac  data, with  good  results. 

2.1.  Introduction 

Detecting  non-rigid  correspondences  between  curves  is  an  important  problem 
with  many  applications  in  computer  vision  and  medical  imaging,  including  motion 
analysis  [56,  15,  45,  41],  shape  analysis[50,  23],  feature-based  registration[51]  and 
knowledge-based  segmentation[72,  13].  Given  a pair  of  curves  Ci  and  C2,  the  problem 
is  formulated  as  finding  a meaningful  transformation  $5  : Ci  — >•  C2  that  minimizes 
a shape-dissimilarity  measure  between  the  curves.  Any  work  on  the  correspondence 
problem  addresses  the  following  issues: 
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2.1.1  Shape  Representation 

Many  models  [58,  83,  81]  represent  the  curves  in  parametric  form  and  find  a 
transformation  between  their  parameter  spaces.  Even  though  these  methods  have  the 
advantage  of  being  computationally  efficient,  it  is  difficult  to  extend  these  methods 
to  match  curves  of  arbitrary  topology. 

A natural  way  to  represent  curves  of  arbitrary  topology  is  to  embed  them  in 
the  plane  implicitly,  as  0-level  sets  of  a 2-D  function.  A lot  of  work  has  gone  into 
developing  algorithms  for  implicitly-defined  curves/surfaces  in  vision  problems[65]. 
In  this  chapter,  we  propose  a framework  to  solve  the  correspondence  problem  for 
implicitly  represented  curves. 

2.1.2  Modelling  the  Transformation 

Recently,  some  techniques  have  been  proposed  to  match  implicitly-defined 
curves[67,  36].  In  general,  such  methods  find  a suitable  transformation  T : 3?2  -*  5i2 
that  minimizes  the  Euclidean  distance  between  the  curves.  These  spatial  transforma- 
tions T can  be  modeled  using  splines[7,  9]  or  dense  deformations[67].  However  in  this 
case,  the  cost  of  accurately  computing  local  variations  between  the  curves  is  high, 
since  it  has  to  be  achieved  by  a transformation  of  the  embedding  domain.  Moreover, 
this  does  not  guarantee  the  point-correspondence  criterion  T(Ci)  = C2,  which  is  es- 
sential in  many  applications,  (e.g., motion-tracking).  These  methods  are  more  suited 
for  curve-matching,  where,  the  goal  is  to  compute  the  shape-dissimilarity  between 
curves  and  hence,  a point-correspondence  is  not  needed.  In  contrast,  our  work  is 
applicable  both  to  curve-matching  and  to  motion-tracking  of  implicitly  represented 
curves. 

One  of  the  novelties  of  our  work  is  that  we  solve  the  correspondence  problem 
for  implicitly  represented  curves  by  directly  computing  transformations  between  the 
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curves.  Here,  worthy  of  mention  is  the  work  of  Memoli  et  al.  [59],  to  find  maps 
between  two  general  curved  manifolds. 

To  be  more  precise,  we  look  for  a diffeomorphism  S : Ci  ->■  C2.  Recently 
there  has  been  interest  in  modeling  diffeomorphic  transformations  for  landmark 
matching[38,  31,  55]  and  image  matching[25]  using  flows  of  time-dependent  vector 
fields.  This  approach  has  been  found  to  be  suitable  for  generating  large  non-rigid 
deformations  of  the  plane.  We  use  a similar  idea  to  construct  our  search  space  of 
diffeomorphisms  between  the  curves. 

To  simplify  computations,  we  map  our  curves  initially  to  a topologically  equiv- 
alent object,  say  C(e.g.,  in  the  1-curve  case,  C is  chosen  to  be  a circle.).  We  then  look 
for  diffeomorphisms  of  C.  A similar  approach  is  taken  in  [31,  2],  but  for  matching 
landmarks  on  surfaces. 

2.1.3  Shape-Matching  Criterion. 

The  similarity  measure  chosen  for  a shape-matching  problem  usually  depends 
on  the  type  of  variability  expected  in  that  class  of  shapes.  Euclidean  distance  [67,  36], 
curvatures [81,  45]  and  curve  normal[41]  are  some  of  the  criteria  that  can  be  used.  For 
a survey  on  similarity  measures,  we  refer  the  reader  to  [86].  The  method  we  introduce 
in  this  chapter  is  geared  towards  curves  that  differ  from  each  other  primarily  by  local 
stretching  and  piecewise-rigid  motions.  For  this  type  of  variability,  we  propose  to  use 
curve  normals  as  the  shape-matching  criterion. 

2.2.  Related  Work 

In  the  work  of  Glaunes  et  al.  [30],  curve-  and  surface-matching  are  treated 
as  specific  cases  of  diffeomorphic  matching  of  measures  on  -ft2  or  SR3.  The  point- 
matching algorithm  of  H.  Guo  et  al.  [32]  can  be  used  for  matching  point  sets  on 
curves  and  surfaces.  These  methods  search  for  spatial  diffeomorphisms  that  minimize 
the  distance  between  the  point  sets  and  have  the  ability  to  robustly  reject  outlier 
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points.  However,  neither  of  these  techniques  uses  shape  as  a matching  criterion, 
which  is  crucial  for  modeling  non-rigid  variations.  Further  we  use  a reduced  search 
space  of  diffeomorphisms  between  the  curves,  rather  than  that  of  the  whole  plane. 
This  eliminates  the  need  for  any  constraint  in  the  objective  functional  to  enforce 
T(Ci)  = C2.  Moreover,  we  can  generalize  our  framework  to  track  curves  of  different 
topology. 

Our  work  can  be  seen  as  a modification  and  extension  of  Tagare’s  shape- 
matching technique  [81]  to  implicitly-defined  curves.  Tagare  et  al.  use  a specific 
class  of  correspondences  called  bimorphisms  to  match  curves.  We  briefly  review  his 
work  here. 

Let  Ci  and  C2  be  two  curves  with  lengths  LCl  and  L c2,  parametrized  by 
arc-length.  A parametrized  bimorphism  is  a map  fi  : [0, 1]  — l [0, 1]  x [0, 1],  given 
by  = tlT7>  r|:]>where  Si(t)  : [0, 1]  -*  [0,LCl]  and  s2(t)  : [0,1]  ->  [0,Lc2]  are 
surjective  functions,  satisfying  the  conditions  s^ (t)  > 0,  s^t)  > 0 and  |/i(t)|  > 0. 
The  functions  Si(t)  and  s2(t)  thus  defined  intrinsically  model  continuous,  non-rigid 
correspondences  between  the  curves. 

Tagare  et  al.  define  a shape-based  correspondence  between  Ci  and  C2  to  be  a 
bimorphism  /i(t)  = [^^,  ^p]  that  minimizes  the  following  energy: 

Jc„o,[s1.s2]=  f [lA(Sl(t))-l®i(s2(t))]2dS„  (2.1) 

J[p]  clSp 

where  0;(s)  is  the  angular  orientation  of  the  normal  to  the  curve  Q at  arc  length  s 

and  dsM  = * + ^r-  dt  is  the  arc-length  element  of  [ptj:=  Image(/x). 

y Cj  c2 

The  objective  function  J can  be  seen  to  be  symmetric  and  scale/rotation- 
invariant  with  respect  to  Ci  and  C2. 

Since  what  we  want  to  minimize  in  (2.1)  is  (^[©i(si(t))  - 02(s2(t))])2,we 
can  use  one  derivative  fewer  in  our  formulation  by  seeking  a constant  6 such  that 
0i(si(t))  — ©2 (s2 (t))  ~ 9.  To  use  a similar  idea  to  match  implicitly-defined  curves,  it 
is  convenient  to  compare  the  curve-normal  functions  ni(si(t))  and  n2(s2(t))  instead. 
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This  motivates  us  to  reformulate  (2.1)  as 


Jcltc2[si,s2,0]  = / MsiOO)  - Ren2(s2(t))|2dsM 

Jm 

where  Rg  is  the  counter-clockwise  rotation  by  the  angle  6. 


Jc,,c2[si»s2,0]  = / 


(2.2) 


2.3.  Bimorphism  for  Level  curves 


Let  and  fl2  be  bounded  regions  in  K2.  Given  Lipschitz  functions 
: fli  — > 5?  and  $2  : ~ > 3?  with  0-level  curves  Ci  and  C2  , we  extend  the  problem 

of  computing  shape-based  bimorphisms  for  these  embedded  curves.  We  require 
to  have  non- vanishing  gradients  on  their  0-level  sets.  E.g.,  <f>j  can  be  chosen  as  the 
signed  distance  function  of  Cj,  in  a small  neighborhood  of  Q and  then  smoothly 
extended  to  the  rest  of  fl;.  For  the  sake  of  simplicity,  we  start  with  the  case  when 
Ci  and  C2  have  the  topology  of  a circle.  Let  C be  a unit  circle  in  -ft2  centered  at  the 
origin,  and  $ : 12  -*  5i  its  signed  distance  function. 


We  define  a bimorphism  between  the  embedded  curves  Ci  and  C2  to  be  a map 
v : C — > (Ci  x C2)  C K4  given  by: 


Lg;  = 5($j)| V$;|dx  is  the  length  of  C;.  Here  <5($(x))  is  the  Dirac  delta  func- 

tion. In  our  implementation,  we  use  a smooth  approximation  for  the  delta  function, 
denoted  5e(t). 

Along  the  same  lines  of  (2.2),  in  this  paper  we  define  a shape-based  correspon- 
dence between  Ci  and  C2  to  be  a bimorphism  that  minimizes  the  energy: 


where  f;  : C — > Q is  a diffeomorphism. 


where  v;  = j^lq,  is  the  unit  normal  vector  field  to  the  level  sets  of  and  ds^,  is  the 


length  element  of  v. 


Let  Ne(C)  = {x  G f2|$(x)  < e},  N£i(Ci)  = {x  £ fli|$j(x)  < ei}  be  neigh- 
borhoods of  C and  C;  respectively,  such  that  we  can  extend  the  diffeomorphism 
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fi,f2 

/ \ 


Figure  2.1:  The  shaded  gray  bands  are  Ne(C),  N£i(Q). 
fo,i  are  initial  diffeomorphisms.  f(  are  the  minimization 
variables. 

fi  : C ->■  C;  to  a diffeomorphism  fi  : Ne(C)  ->  N€.(Cj),  i.e.  fi|c  = Consider  the  map 
F : Nc(C)  — >■  5ft4  defined  as  F = [j^-,  j^-].  The  graph  of  the  map  F is  a surface  in 
5ft4,  with  surface-element  |FX  A Fy|dx  = A/|Fx|2|Fy|2  - |FX  • Fy|2dx. 

Since  F|c  = v,  the  length  of  the  curve  v embedded  in  5ft4  can  be  computed  by: 

[dsv=  [ 5(<f>)|V<f>||Fx  A Fy|dx  = [ <5($)|FxAFy|dx, 

Jv  J n Jq 

since  |V$|  = 1.  Hence,  we  can  rewrite  (2.3)  in  the  form 

E(fi,f2,*)=  f ^(^(x))|v!(fi(x))  - R0v2(f2(x))|2|Fx  AFy|dx  (2.4) 

J n 

2.4.  Decomposition  of  the  Bimorphism 

In  the  energy  (2.4),  the  search-space  is  over  diffeomorphisms  fi,f2  : Ne(C)  — > 
5ft2, with  fi(C)  = Q. 

Let  f0)i  : Ne(C)  — > f2i(i=l,2),  be  diffeomorphisms  arbitrarily  chosen  such  that 
fo,i(C)  = C;.  Then  given  fj  there  exists  a unique  diffeomorphism  f;  : Ne(C)  — > Ne(C) 
with  fj(C)  = C,  such  that  f;  = f0!;  o f;(Fig.2.1). 

Thus  we  can  replace  our  search-space  by  the  space  of  diffeomorphisms  fl5f2  : 
Ne(C)  — > Ne(C)  such  that  fi(C)  = C.  Now  the  energy  functional  (2.4)  is  written  as 
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E(fi,f 2,0)=  f <K$(x))|vi(f0>i  ofi)  - Rev2(fo,2  of2)|2|Fx  AFy|dx  (2.5) 
Jn 

where  F = 

We  can  observe  that  the  set  of  minimizers  to  energy  functional  (2.5),  if  non- 
empty, is  invariant  with  respect  to: 

1.  Initial  maps  f0,j. 

2.  Rotation  and  scaling  transformations  on  the  level  curves  Ci  and  C2. 

To  search  for  f;,  we  consider  vector  fields,  only  of  the  form  X,  = aj(x)T,  where 
T = V^-1-  is  the  tangential  vector  field  to  the  level-curves  of  $ obtained  by  counter- 
clockwise rotation  of  V<3>,  and  where  a\  : f)  — » 5R  is  a Lipschitz  function  with  compact 
support  in  Ne(C). 

For  a fixed  x,  Consider  the  flow  Xi(x,  t)  of  X;  given  by  the  following  differential 
equation: 

Xi(x,t)  = Xi(xi(x,t)), 

Xi(x,0)=x.  (2.6) 

Then  for  a fixed  r > 0,  Xi(xi T)  '■  Ne(C)  — > N£(C)  is  a diffeomorphism  and  Xi(C,  t)  = 

C.  Naturally,  we  let  fi(x)  = fi(x;  a;)  = Xi(x,  T)-  Hence  our  search  for  diffeomorphisms 

fi  : Ne(C)  — > Ne(C)  satisfying  fi(C)  = C,  reduces  to  looking  for  Lipschitz  functions  a; 

with  compact  support  in  Ne(C).  The  energy  functional  (2.5)  now  becomes: 

E0(ai,  a2,  9)  = f <5($(x))|vi(f0)i  o fx(x;  ax))  - R0v2(fO)2  o f2(x;  a2))|2 
J n 

|FX  A Fy|dx  (2.7) 


2.5.  Approximation 

We  approximate  our  minimization  problem  (2.7)  by  restricting  our  search 
space  for  ai  and  a2  to  be  finite-dimensional.  Let  {xj}]!!  be  a finite  set  of  regularly 
spaced  points  on  C.  We  look  at  N basis-functions  with  ipj  centered  at  the 
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control-points  xj,  defined  by  ^j(x)  = ?7£(r)(l  - cos  1(cj))  , where  r = |x|  - l,q  = 

7T  J I X I 

rje  £ C£°(5i)  is  the  standard  mollifier  of  width  e. 

Let  d'(x)  = [^i(x),  ^(x), ^n(x)].  To  solve  our  problem  , we  find  constants 
ci  = [cii,ci2,  ...,c1N]  and  c2  = [c2i,  c22,  • c2n],  and  such  that  a;  = = Cj-'L, 

and  6 minimize  (2.7)  over  !ftN  x JiN  x 3?/27tZ.  Hence  a;  £ C°(N£(C)),  as  needed. 

Since  the  energy  E0  depends  only  on  the  restriction  a;|c,  a;  can  be  computed 
in  any  manner  away  from  C,  so  that  a;  £ C°(Nt(C)).  Hence,  all  that  we  need  to 
consider  while  choosing  is  that,  any  element  of  C°(C)  can  be  approximated  to 
desired  accuracy,  by  linear  combinations  of  {^j|c}jli-  For  our  choice  of  fa,  it  is  seen 
that  the  approximating  functions  a,  are  piecewise  linear  on  C and  for  a sufficient 
number  of  control  points,  aj|c  can  closely  approximate  any  member  of  C°(C).  The 
mollifier  rje  ensures  an  annular  support  of  width  e for  a;. 

To  summarize,  we  minimize: 

EN(ci, c2, 6)  = [ 5($(x))|v1(f0,i  of^xsci))  - Rev2(fo,2  of2(x;c2))|2 
Jn 

|FxAFy|dx  (2.8) 

over  3?N  x»Nx  3i/27rZ. 

Here  F = fj(x)  = Xi(x,  r)  where  Xi  satisfy  the  flow  equations: 

Xi(x,t)  = X;(xi(x,  t)), 

Xi(x,  0)  = x.  (2.9) 

For  any  map  T,  we  denote  its  Jacobian  matrix  as  DT.  Denote  fic.,  = Let 
V = Vj (fo,i  o fx(x;  Ci))  - Rev2(f0)2  o f2(x;  c2)).  The  Euler-Lagrange  equations  for  this 
minimization  problem  are: 

£En  r f) 

fr-  = / Wx))V.Dv1.Df„,1.f1,Clk)|R,  A Fy|  + |V|2—  (|FX  A Fyj)dx  = 0,  (2.10) 
c'cik  Jn  o Cik 

0£N  r a 

7j—=  / i(l>(x))(-V.Ri,Dv2.Df0i2.f2,ejk|F,AFJ|  + |V|2— (IF.AFjDdx 

&c2k  Jn  00,2k. 

— 0,  (2.11) 
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(2.12) 


9EN  r 

-gf=  [f«x)){-V.Riv!)!F.flF,|=0. 

In  (2.10)  and  (2.11), 

a^(lFx  AFy|)  = Tf^F^f[|Fy|2Fx-FXCil(  + |F*|2Fy -FyCjk  - (Fx-Fy)(Fx-FyCik  + Fy  *Fxc.fc)]. 
To  simplify  further,  we  consider  fijCik(x; Ci).  Since  f;(x)  = Xi(x,r),  we  can  compute 
fi,Clk  (x)  = Xi,Cik(x,  r)  from  the  differential  equation: 

Xi,'cik(x,  t)  = V’k(Xi)T(xi)  + (Ci  • T(xi))DT(xi)  • Xi,cik, 

Xi,cik  (x,  0)  = 0. 

Using  this  in  F(x;c„c2)  = [Su^SSl,  Sfigial],  we  get 


(2.13) 


Clk 


Dfoa  ' fl’Clk , 0]  and  FC2k  = [0,  ?f°’2  ' f;2’C2k 


jc2 


]• 


Since  FXCik  = FCikX  and  Fyc.k  = FCiky,  we  can  write: 

af^(lFXAFy|)  = ]F^FU[|Fy|2FX'FCikX+|Fx|2Fy  -F^yy-  (Fx-FyXFx'Fc^y  + Fy-Fqyy)]. 


2.6.  Numerical  Implementation 

To  discretize  Eqs.(2.10),  we  use  a finite  difference  scheme.  For  any  map  T 
defined  on  let  Tpq  denote  its  value  at  the  pixel  xpq  = [p/j,  qh],  where  h is  the  pixel 
size.  We  use  gradient-descent  to  minimize  the  energy  (2.8)  over  5RN  x 5?N  x $R/27rZ. 
Given  initial  values  cj,  c°  and  6°,  we  iteratively  solve  for  c™,  c™  and  9m,  (m=l,2,3....). 

For  an  iterate  cf1,  9m,  let  us  denote  f^pq  = fi(xpq)  = Xi(xpq,r).  We  can 
compute  f^q  by  numerically  solving  the  equation: 

Xi(xPq,  t)  = (c,m  • T(xi(xpq,t)))T(xi(xpq,t)), 

Xi(Xpq,0)  = Xpq.  (2.14) 

for  t e [0,  t], 

Let  fi^ik,pq  denote  ^(Xpqwr1)  = £(xpqir)>  (k  = 1,  --N).  We  can  compute 
f™ik)Pq  by  numerically  solving  the  equation: 
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(2.15) 


for  t £ [0,  r].  At  the  iteration  m,  the  bimorphism  is  given  by  Fm  = [f°^°fl  , f°L2°f^  ]. 
Hence,  writing  um  = [c™,  c™,  9m]  we  arrive  at  the  discrete  form  for  the  Euler-Lagrange 
equations  (2.10)-(2.12): 


(k  = 1,  ..N),  where  Je(t)  is  a regularized  approximation  for  <5(t)(see  [9]  for  details).  In 
(2. 14)-(2. 18) , the  derivatives  are  approximated  using  first-order  difference  schemes. 
We  summarize  the  algorithm  as  follows: 

Input  initial  u° 
for  N = 1,4,8... 
for  m = 0, 1, 2... 

1.  From  um  obtain  (fim,f^,f^lk,f£cJ  using  (2.14,2.15). 

2.  Compute  VEN(um)  using  (2.16)-(2.18). 

3.  Update  um+1  using  gradient  descent. 

end 

Re-initialize  u°. 

end 


p>q 


A F"|)  = 0,  (2.16) 


(2.18) 


(2.19) 
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Figure  2.2:  f0ji  is  computed  by  mapping  the  discrete  set 
of  points  on  N£(C)  to  those  in  the  neighborhood  of  Q 
and  then  interpolating  to  the  whole  of  Ne(C). 

The  step  N=l,4,8...  is  necessary  for  the  algorithm  to  avoid  local  minimizers 
of  (2.7).  We  start  with  a given  u°  for  the  N=1  step  and  re-initialize  u°  for  next  value 
of  N,  by  interpolating  the  minimizer  for  the  current  value  of  N.  This  is  continued 
until  the  change  in  energy  (2.7)  goes  below  a threshold.  The  diffeomorphisms  f0,i 
and  f0j2  map  the  annular  region  Ne(C)  to  neighborhoods  of  Ci  and  C2(Fig.2.1).  Note 
that  the  minimum  value  (if  it  exists)  is  independent  of  choices  of  f0)i.  Hence  these 
maps  should  be  computed  initially  in  the  simplest  possible  manner,  such  as  by  the 
following  procedure: 

Let  {Pj)}jl1  and  {aj)}jl1  be  evenly  spaced  points  on  the  curve  Ci  and  the  circle 
C,  placed  counter-clockwise  along  Ci  and  C with  increasing  j.  Let  fo,i(aj))  = pjh  We 
now  extend  f0,i  to  the  neighborhood  Ne(C).  For  j=l,2...M,  let 
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n 


Figure  2.3:  Directed  Graphs  for  level  curves.  II  The 

dots  represent  the  vertices  and  the  directed  line  seg- 
ments represent  the  directed  edges. 


Pj  1 = Pj°  - (1  + £*j)nj, 

P}  = P?  + (!  + (2.20) 

where  acj  is  the  curvature  and  nj  = is  the  unit-normal  vector  of  the  curve  Ci  at 
the  point  p°.  £ > 0 is  a constant  to  be  chosen  below. 

Similarly,  for  the  circle  C,  let 

af1  = a?  - eilv 

aj  = a?  + enj,  (2.21) 

where  nj  = is  the  unit-normal  vector  of  the  circle  C at  the  point  a?. 

For  k = -1,1  and  i = 1,2,....M,  we  define  fo^aj5)  = pjb  (Fig.2.2).  We  can  now 
interpolate  f0,i  with  as  much  smoothness  as  required  to  other  points  in  Nf(C).  We 
take  £ sufficiently  small  that  f0  is  a diffeomorphism.  The  initial  diffeomorphism  f0j2 
is  computed  in  a similar  way. 
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2.7.  Curves  of  Arbitrary  Topology 

So  far  we  have  restricted  our  attention  to  the  case  when  the  curves  being 
matched  have  the  circle-topology.  This  model  can  be  used  to  match  two  arbitrary, 
closed,  0-level  curves  of  the  same  topology.  An  interesting  extension  is  to  match 
closed,  O-level-curves  with  different  topologies.  Below,  we  briefly  introduce  shape- 
based  bimorphisms  between  implicitly-defined  curves  Ci  and  C2  of  arbitrary  topology. 

2.7.1  Topological  Configuration  of  Level  Curves 

In  this  thesis,  we  assume  that  the  level  curves  we  are  dealing  with  are  dis- 
joint unions  of  simple  closed  curves.  Given  such  a level  curve,  we  can  describe  its 
topological  configuration  by  means  of  a directed  graph.  A directed  graph  is  a pair 
< V,  E >,  where  V is  a set  of  vertices  and  E C V x V is  a set  of  directed  edges. 
For  any  (r,  s)  G E,  r is  the  initial  vertex  and  s,  the  end  vertex.  So,  to  a level 
curve  C = Ul=i  where  each  is  a single  curve,  we  can  associate  a directed 
graph  <V,E>,  where,  V = {1,2,  ...,T}  and  E consists  of  all  pairs  ( k,l ) such  that 
Int(C^) nExt(C^)  — C is  connected. (Fig. 2. 3)  shows  a level  curve  with  its  associated 
directed  graph. 

2.7.2  Multiple-Curve  Case  of  Same  Topology 

Let  us  suppose  that  Ci  and  C2  have  the  topology  of  the  union  of  T disjoint 
circles,  both  with  the  same  directed  graph  G.  Let  C = (J^=1  C(a),  each  C(q)  a circle, 
have  a topology  consistent  with  G,(Fig.2.4).  If  we  choose  e so  that  the  neighborhoods 
N e(C<“>)  are  mutually  disjoint,  then  we  can  find  smaller  neighborhoods  Nei(Ci)  of 
Cj,  so  that  any  diffeomorphism  : C — > Q can  be  extended  to  a diffeomorphism 
f;  : Ne(C)  -4  Nei(Q).  Given  initial  diffeomorphisms  f0,i  : Ne(C)  -4  Nei(Q),  we 
look  for  Lipschitz  functions  ai,a2  and  a locally  constant  function  0,  i.e  one  which  is 
constant  on  each  connected  component  of  Ne(C),  that  minimize  an  energy  functional 
similar  to  (2.7): 
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Figure  2.4:  Framework  for  computing  bimorphisms  for 
multiple-curves  of  same  topology. 


Ei(aj,  a2,  0,  a)  — f 5($(x))|v(f0  o f(x;  a))  - R0w(go  o g(x;  b))|2 
Jn 

|FxAFy|dx  (2.22) 

The  function  0 recovers  relative  rotations  between  each  of  the  T corresponding 
components  of  Ci  and  C2.  o runs  over  the  symmetries  of  G,  the  directed  graph  of 
C(e.g.,  the  symmetry  group  of  the  graph  in  Fig.2.3  is  Z2  x Z2). 

Minimizing  (2.22)  is  the  same  as  matching  the  components  of  Ci  and  C2  one  at 
a time,  using  (2.7).  However,  one  can  easily  add  spatial  constraints  for  the  functions 
ai,a2  and  0 within  (2.22),  thus  cohering  the  deformations  of  the  components  of  Cj. 
For  this,  a,  and  0 need  to  be  extended  from  Ne(C)  to  Q (which  can  be  done,  for 
instance,  by  approximating  aj  and  0 to  be  linear  combinations  of  spatially  varying 
basis  functions,  as  shown  in  section  2.5). 
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Figure  2.5:  Computing  bimorphisms  between  2 

implicitly-defined  curves  of  different  topology.  The 
points  p1;  p2,  q:  and  q2  can  be  chosen  arbitrarily. 


2.7.3  Multiple-Curve  Case  of  Different  Topology 

For  simplicity,  in  this  thesis  we  consider  only  the  case  in  which  Ci  has  the 
circle-topology  and  C2  has  the  topology  of  two  non-nested  circles. 

Let  qx  and  q2  be  two  points,  one  on  each  of  the  connected  components  of 
C2.  We  denote  C2  = C2  — {q1,q2}.  Let  C be  a circle  and  let  p4  and  p2  be  any 
two  distinct  points  on  C.  We  denote  C = C — {p1?p2}.  A bimorphism  between 
the  implicitly-defined  curves  Ci  and  C2  is  a map  v : C ->  (Ci  x C2)  C given  by 
v = -^r],  where  f2  : C -»  C2  is  a diffeomorphism  and  f:  is  the  restriction  to  C 

of  a diffeomorphism  from  C to  Ci  . (Fig.2.5) 

A shape-based  correspondence  between  Ci  and  C2  is  a bimorphism  that  min- 
imizes the  energy: 


(2.23) 
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Figure  2.6:  Framework  for  computing  bimorphisms  for 
multiple-curves  of  different  topology. 


where  0 is  a locally  constant  function.  It  can  be  seen  that  the  set  of  minimizing 
bimorphisms  for  (2.23)  if  nonempty,  is  invariant  with  respect  to  the  choice  of  pl5  p2, 
Qi  and  q2. 

We  now  construct  a search  space  for  the  diffeomorphisms  fj.  Each  curve  Q 
is  first  mapped  to  a topological  representative  C^(a  disjoint  union  of  circles).  In 
the  previous  sections,  we  saw  a framework  which  can  be  used  to  construct  a space 
of  diffeomorphisms  fj  : -4  Q.  This  framework  was  formed  by  decomposing  f;  as 

fi  = f0,i  o fi;  where  f0,i  : -4  C;  are  initial  maps  and  f;  : C®  -4  C®  are  constructed 

as  time-r  flows  of  tangential  vector- fields  X;  = a;Tj  on  C^.  The  functions  aj,  chosen 
to  be  at  least  Lipschitz,  then  model  the  search  space  for  f. 

Now,  to  search  for  diffeomorphisms  fj,  we  write  f;  = fj  o u;,  where  u;  : C — > 
is  any  diffeomorphism.  The  locally  constant  function  0 used  for  rotations,  is 
defined  over  C(2\  Fig. 2. 6.  Similar  to  previous  cases,  let  Nf(C),  N£(o(C^),  N£.(Q)  be 
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Figure  2.7:  Shape-based  bimorphism  between  2 shark 

curves.  Colored  segments  on  C are  mapped  to  the  same- 
colored  segments  on  C\  and  C2. 


the  corresponding  neighborhoods  of  C,  C(l)  and  Cj  respectively.  Let  fleo  = fl  — 

(Neo(p1)  UNeo(p2)),  where  Neo(pj)  is  a small  neighborhood  of  p;.  The  energy  (2.23) 

can  be  written  as  a minimization  problem  over  the  space  of  Lipschitz  functions  a; 

defined  over  Ne(i)(C^)  and  locally  constant  0 defined  over  Ne(2>  (C^): 

E2(ai,  a2,  0)  = [ 5($(x))|vx(fo,i  o fi(ux(x);  a))  - Rev2(f0,2  o f2(u2(x);  b))|2 
JricQ 

|FX  A Fy|dx  (2.24) 


where,  F = 


fo,l°fl°Ul  fo,2°f2°U2'] 

LCi  Lc2  •* 


2.8.  Matching  Genus-0  Surfaces 

Part  of  our  work  in  progress,  is  to  extend  this  model  to  match  shapes  of  0- 
level  surfaces.  The  basic  framework  remains  the  same,  i.e.surface  normals  are  used 
as  the  shape  matching  criterion  and  diffeomorphisms  are  generated  by  flows  of  vector 
fields.  However,  unlike  in  the  2D  case,  the  sphere  does  not  have  a non- vanishing, 
smooth,  tangential  vector  field.  Hence  we  look  at  flows  of  arbitrary  smooth  vector 
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Figure  2.8:  I.  Initial  correspondence:  colored  segments 
not  matched  properly.  II.  Shape-based  bimorphism 
computed:  Colored  segments  in  Vase  1 and  Vase  2 are 
properly  mapped. 

fields  defined  on  the  sphere,  constrained  to  have  zero  normal  component.  These  flows 
generate  a dense  set  of  diffeomorphisms  of  the  sphere. 

2.9.  Results  and  Applications 

In  this  section,  we  will  demonstrate  the  applicability  of  our  model  to  synthetic 
(Fig. 2. 7,  Fig. 2. 8,  Fig.2.9,  Fig.2.12)  and  cardiac  data  (Fig.2.10).  Given  the  signed 
distance  functions  for  the  two  curves  Ci  and  C2,  we  use  algorithm  (2.19)  to  compute 
a shape-based  bimorphism  v : C — ► Ci  x C2,  v = [f  1 , ^2]  between  the  curves  (see 
Eqn.(2.3)).  The  energy  of  v,  ECl,c2[fi,  4,  0]  gives  a symmetric  distance  between  the 
curves  Ci  and  C2.  Curves  that  differ  only  by  local  tangential  stretching  and  rigid 
motion  are  equivalent  under  this  distance. 

In  (2.7),  we  simplified  Ec1,c2[fii  0]  to  the  form  of  E0[ai,  a2, 0].  The  problem 
is  then  reduced  to  the  minimization  of  EN[c,  d,  0]  = Eo[ai,a2,0],  by  writing  a,  as  a 
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linear  combination  of  N basis  functions  ipj,  with  weights  c and  d.  Algorithm  (2.19) 
minimizes  a sequence  of  functions  EN,  (N  = 1,  2,4..).  The  resulting  correspondences 
are  automatically  diffeomorphisms,  since  the  necessary  smoothness  of  the  flow  ve- 
locities a*,  follows  from  the  choice  of  ipj.  Further,  the  energy  EN  is  symmetric  with 
respect  to  the  order  of  Ci  and  C2-  Since  we  do  not  know  if  the  sequence  of  minimizers 
of  Ex,  actually  converges  to  a minimizer  of  E0,  we  use  experiments(Fig.2.7,  Fig. 2. 8, 
Fig-2.9)  to  validate  the  algorithm  (2.19)’s  ability  to  capture  the  ‘right’  diffeomor- 
phism  between  curves  differing  mainly  by  local  tangential  stretching(possibly  large) 
and  rigid  motions.  Our  experiments  indicate  that  the  sequence  of  minimizers  of  EN 
converge,  and  the  flow  velocities  corresponding  to  this  sequence  decrease  the  energy 
of  the  functional  E0. 

In  our  results,  to  visualize  a bimorphism  between  curves  Ci  and  C2,  the  canon- 
ical domain  C is  partitioned  into  segments,  each  of  a different  color.  Then  each  seg- 
ment’s image  under  the  bimorphism  is  indicated  by  segments  of  the  same  color,  on 
the  curves.  For  instance,  Fig. 2. 7 shows  a shape  based  bimorphism  that  was  computed 
for  two  synthetic  ‘shark’  curves. 

In  Fig.2.8,  we  see  two  ‘vase’  curves  being  matched.  These  curves  differ  only  by 
small,  local  tangential  stretching  and  the  resulting  bimorphism  has  matched  similar 
shapes,  as  expected. 

In  Fig.2.9,  we  test  the  model’s  capability  to  recover  large  non-rigid  deforma- 
tion. We  compute  a shape-based  bimorphism  between  curves  Ci  and  C2  that  differ 
by  a large  tangential  stretching.  Starting  with  an  initial  bimorphism(arbitrarily  cho- 
sen)^), (2.19)  computed  bimorphisms  for  N=l, 4, 6, 8, 10, 12, 14  control  points,  the  final 
result  for  N— 14  shown  in  II. 

We  define  two  error  measures  to  estimate  the  accuracy  of  our  correspondences. 
The  landmark  error (£),  measures  the  discrepancy  between  the  computed  correspon- 
dence and  sets  of  manually  selected  landmarks.  We  chose  landmarks  that  could 
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be  reliably  identified  due  to  their  salient  nature(e.g.,  the  curvature  critical  points). 
Given  L pairs  of  landmark  points  {ii,i2},  C is  defined  as  the  mean  of  the  errors 
Ang[f2°fi_1(ii)  — c,  i2~ c],  where  c is  the  centroid  of  the  curve  C2  and  Ang[x,y]  £ [0, 1] 
is  the  angle  between  vectors  x and  y normalized  by  i r.  Hence  the  maximal  error  of 
100%  occurs  when  the  mapped  point  and  the  landmark  are  180°  apart. 

The  normal  error(n)  measures  the  misalignment  between  the  normals  at  the 
matched  points  on  the  curves.  Given  a set  of  R evenly  distributed  points  pj  on  C,  n 
is  defined  as  the  mean  of  the  errors  ^(^(p;))  - R0n2(f2(pi)|2.  Ideally  for  curves  that 
differ  only  by  local  tangential  stretching  and  rigid  motion,  for  some  fx,f2  and  #,  n=0 
for  all  R. 

Table  2.1  shows  the  errors  ( and  n for  the  stretch  test(Fig.2.9),  during  dif- 
ferent stages  of  the  minimization.  The  energy(E0)  and  the  errors(£  and  n)  are  seen 
to  initially  decrease  rapidly,  and  then  stabilize  for  the  optimal  number  of  control 
points(here  N=14). 


Table  2.1:  Errors  C and  n for  the  stretch  test  for  N=l, 4, 6, 8, 10, 12, 14  control  points. 


N 

Iteration 

Energy  (E0) 

Landmark 

error(£%) 

Normal 
error  (II) 

1 

5 

1639 

13.09 

1.37 

4 

9 

923.07 

10.13 

.59 

6 

30 

459.45 

4.05 

.230 

8 

51 

374.3 

4.84 

.239 

10 

72 

336.78 

4.28 

.214 

12 

93 

312.54 

3.78 

.194 

14 

114 

300.54 

3.85 

.201 

Since  local  contraction/stretching  is  the  main  component  in  cardiac  motion, 
curve-normals  is  a viable  matching  criterion.  In  Fig.2.10,  using  Tagare’s  data  (avail- 
able in  his  website  http://noodle.med.yale.edu  / hdtag /)  of  the  endocardium  of  a 
dog’s  heart,  we  find  a shape-based  bimorphism  between  the  boundaries  of  the  heart 
at  end-systole  and  end-diastole.  In  (I-III) , we  see  the  correspondence  during  the 
stages,  N=4,8,16  of  algorithm  (2.19).  We  picked  pairs  of  key  landmark  points  on 
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Table  2.2:  Errors  C and  II  for  the  test  on  cardiac  data  for  N=l,4,8,12,16  control 
points. 


N 

Iteration 

Energy(Eo) 

((%) 

n 

Ansi 

ANa2 

1 

6 

395.45 

19.02 

.82 

- 

- 

4 

17 

176.24 

6.22 

.39 

4.42 

.86 

8 

47 

137.15 

4.01 

.36 

.21 

1.80 

12 

57 

120.82 

3.55 

.31 

.88 

.50 

16 

74 

120.62 

3.55 

.31 

.0018 

.0032 

the  outer(circles)  and  inner(boxes)  curves.  This  was  used  as  the  ground  truth  for 
the  landmark  error (C).  The  black  lines  connecting  the  points  on  the  outer  and  inner 
curves,  denote  the  correspondence  that  was  obtained  by  our  algorithm  for  the  same 
landmark  points  (circles)  on  the  outer  curve. 

To  minimize  (2.7),  we  start  with  N=1  in  (2.19)  with  c=O,d=O,0  = 0,f  = 
.4,  ci  = 1.  We  found  stable  convergence  for  N=16  in  74  iterations.  The  graph  in 
Fig.2.11  plots  the  energy  E0(ai,a2,$)  versus  the  number  of  iterations.  The  black 
squares  show  the  minimum  that  was  obtained  for  each  N.  The  algorithm  stops  when 
the  change  in  the  ads  for  consecutive  N (ANa;)  falls  below  a threshold. 

Table. 2. 2 shows  the  errors  ( and  II,  along  with  the  quantities  A^a;,  during 
different  stages  of  the  minimization.  The  quantities  A>jaj  rapidly  decreases  towards 
zero,  and  the  energy  E0  along  with  the  accuracy  errors  (£  and  II)  decrease  initially 
and  stabilize  at  N=16.  This  opens  up  the  possibility  of  proving  the  convergence  of 
the  minimizers  of  EN,  to  a minimizer  of  E0.  This  is  a work  in  progress. 

In  Fig.2.12,  we  used  the  models  presented  in  (2.22)  and  (2.24)  for  match- 
ing/tracking curves  of  arbitrary  topology,  is  shown.  The  results  shown  are  prelimi- 
nary and  are  on  synthetic  data.  In  (I),  the  curve  shown  on  the  left  splits  into  two 
curves.  Colored  segments  of  similar  shapes  are  tracked  before  and  after  the  splitting, 
using  (2.22).  (II)  and  (III)  show  matching  of  multiple  curves  of  the  same  topology. 
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Table  2.3:  Comparisons  of  our  method  with  several  current  implicit  matching  meth- 
ods 


Methods/ 

Features 

Our  Method 

Paragios  [67] 

Guo  [32] 

Glaunes[30] 

Applications 

Motion  Track- 
ing, Shape 

Matching 

Shape  Match- 
ing 

Shape  Match- 
ing 

Shape  Match- 
ing 

Outliers 

Rejected? 

No 

Yes 

Yes 

Yes 

Shape  Based? 

Yes 

No 

No 

No 

Large  de- 

formations 
allowed? 

Yes 

No 

Yes 

Yes 

Search  Space 

Diffeos  between 
curves 

Smooth  trans. 
between  2D  do- 
mains 

Diffeos  between 
2D  domains 

Diffeos  between 
2D  domains 

Symmetry  in 

C1&C2  ? 

Yes 

No 

No 

No 

Invariance 

Scale/Rotation 

Scale/Rotation 

Rotation 

Rotation 

Topology 

Arbitrary 

Same  for 

C!&C2 

Same  for 

C1&C2 

Same  for 

Ci&C2 

2.10.  Comparison 

Table  2.3  shows  a comparison  of  the  model’s  capabilities  with  some  of  the 
existing  implicit  matching  methods. 

Unlike  existing  methods  that  are  suitable  for  curve-matching,  our  model  is 
applicable  to  both  shape-matching  and  to  motion-tracking.  By  definition,  the  prob- 
lem is  to  find  diffeomorphic  correspondences  between  given  two  0-level  curves.  The 
input  curves  Ci  and  C2  are  required  to  be  at  least  Lipschitz  continuous  and  we  do 
not  handle  the  case  of  data  with  outliers. 

Except  for  Paragios’  method,  the  other  methods  propose  a large  diffeomor- 
phism  setting  using  flows  of  suitable  vector-fields  on  the  entire  image  domain.  How- 
ever we  search  for  diffeomorphisms  just  between  the  curves,  thus  resulting  in  a smaller 
search  space  for  the  vector  fields. 

Further,  our  use  of  a shape-based  matching  criterion(i.e. normals)  makes  the 
energy  invariant  under  scaling  and  rotation  of  the  curves  Ci  and  C2. 
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Also  our  energy  functional  is  minimized  over  a “lifted  space”  of  bimorphisms, 
a set  of  curves  on  a torus  Ci  x C2,  satisfying  certain  conditions.  This  approach  makes 
the  resulting  correspondences  symmetric  under  interchange  of  Ci  and  C2. 

Lastly,  unlike  the  mentioned  methods,  we  have  a natural  way  of  match- 
ing/tracking curves  of  arbitrary  topology. 

2.11.  Summary  and  Conclusions 

Given  two  signed  distance  functions,  we  have  presented  a variational  approach 
to  compare  the  shapes  of  their  0-level  curves,  which  has  applications  in  both  shape- 
matching and  motion-tracking  of  implicitly  represented  curves.  With  an  assumption 
that  the  curves  of  interest  differ  primarily  by  piecewise  rigid  motions  and  local  stretch- 
ing, curve-normals  were  used  as  the  shape-matching  criterion.  We  gave  a precise 
definition  for  a shape-based  correspondence  between  the  implicitly-defined  curves, 
specifically  a bimorphism  that  minimizes  the  shape-matching  measure  (Eqn.3).  Even 
though  most  of  the  presentation  focused  on  the  problem  of  matching  two  curves  of 
circle-topology,  the  model  admits  a natural  extension  for  constructing  shape-based 
bimorphisms  between  curves  of  arbitrary  topology.  Further  the  extention  of  the 
framework  to  3 dimensions,  as  described  in  section  8,  looks  promising. 
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Figure  2.9:  Large- deformation  test.  I.  Initial  correspondence.  II.  Shape-based 
bimorphism  computed. 
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Figure  2.10:  Tagare’s  MRI  Heart  Data  I II  & III.  correspondence  for 

N=4,8,16  control  points. 


Figure  2.11:  Convergence  history  for  heart  data(Fig.2.10).  The  squares  show 
the  minimum  that  was  obtained  for  (N=l,4,8,12,16)  control  points. 
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Figure  2.12:  Matching  level  curves  of  arbitrary  topology 
I.  Tracking  splitting  curves.  II  and  III.  Multiple  Curve 
matching. 


CHAPTER  3 

PDE  BASED  METHOD  FOR  FUZZY  SEGMENTATION 


We  propose  a novel  variational  approach  to  automatically  segment  medical 
images  into  a fixed  number  of  classes.  These  classes  form  a partition  of  the  image 
domain  into  connected  regions.  They  are  defined  such  that  intensity  values  within  a 
class  are  close  to  a constant. 

Clustering  methods  such  as  the  fuzzy  c-means(FCM)  are  known  for  their  com- 
putational speed  and  ease  of  implementation.  Recently,  active  contours  have  become 
popular  due  to  their  consistent  mathematical  formulation  and  robustness  to  noise, 
spurious  edges  and  boundary  gaps. 

Our  method  enables  us  to  use  fuzzy  classification  (defined  later)  and  active 
contours  in  a single  variational  framework.  This  approach  allows  the  use  of  tools 
from  both  deformable  geometry  and  clustering  in  a well-defined  setting  and  provides 
a useful,  unsupervised  segmentation  technique.  Further,  this  framework  can  be  ex- 
tended to  segmentation  of  local  structures. 

The  model  was  tested  on  synthetic  and  MRI  brain  data,  with  promising  re- 
sults. 

3.1.  Introduction 

Image  segmentation  plays  a crucial  role  in  many  medical  imaging  applications 
such  as  visualization  of  anatomical  structures,  detection  of  pathology  and  computer- 
assisted  surgery.  The  problem  is  to  partition  the  image  domain  into  disjoint,  con- 
nected components  that  are  homogeneous  with  respect  to  some  image  feature  (e.g., 
intensity,  texture). 


35 


36 


That  is,  given  an  image  I with  intensity  as  a feature,  defined  on  a rectangular, 
closed  region  fi,  we  want  to  find  mutually  disjoint  regions  f 2k  C Q,  k=l,2...C  such  that 
Uk=i  ^ and  each  region  {f2k}k=i,2,..,c  is  homogeneous  with  respect  to  image- 
intensity.  A related  problem  is  image  classification  which  deals  with,  choosing  a 
unique  label  for  each  pixel  from  a fixed  set  of  labels  or  classes  that  are  characterized  by 
image  intensity.  Partitioning  the  image  domain  into  sets,  whose  pixels  belong  to  the 
same  class  is  segmentation.  For  most  of  the  medical  imaging  modalities,  segmentation 
is  a hard  problem  because  of  the  presence  of  noise,  artifacts,  low-contrast  and/or 
poor-resolution. 

In  medical  imaging  modalities  such  as  MRI  and  CT,  the  sampling  process 
results  in  partial-volume  effects  leading  to  structural  ambiguities.  Partial-volume 
effects  are  artifacts  that  occur  when  multiple  tissue  types  contribute  to  a single  pixel, 
resulting  in  a blurring  of  intensity  across  boundaries.  A common  way  to  handle  this 
problem  is  to  allow  for  the  segmented  regions  or  classes  to  overlap,  termed  as  soft 
segmentation.  In  ‘Hard’  segmentation  techniques,  a pixel  either  belongs  to  a class 
or  not.  In  soft  segmentation,  a pixel  is  allowed  membership  degrees  with  respect  to 
each  class,  which  allows  for  uncertainty  in  the  segmentation  process  thus  retaining 
more  information  than  hard  segmentation  methods. 

Hard  segmentation  amounts  to  finding  characteristic  functions  Xk,  correspond- 
ing to  the  regions  k=l,2,..,C.  In  soft  segmentation,  we  look  for  C smooth  functions 
Uk  : — > [0,1],  the  membership  function  for  the  class  k,  (k=l,2..,C).  Fuzzy  tech- 

niques are  soft  segmentation  methods  that  impose  the  partition  of  unity  constraint 
on  {uk}k=i)2..c; 

c 

X>  = 1 

k=l 

The  functions  {uk}k=i,2..c  are  soft  analogies  to  the  characteristic  functions 
{Xk}k=i,2..c-  Just  as  the  value  Xk(x)  indicates  if  the  pixel  x G or  not,  uk(x)  gives 
the  probability  of  membership  of  x in  the  class  k.  If  required,  a hard  segmentation 
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could  be  computed  from  membership  functions  such  as  the  maximum  membership 
value,  i.e  Dk  = (x|uk(x)  > Ui(x)|¥k)1<j<c}. 

In  our  work,  we  look  at  a fuzzy  segmentation  approach  posed  in  a variational 
framework,  with  image-intensity  as  the  distinguishing  feature.  However,  it  is  possible 
to  extend  the  model  to  handle  multi-channel  data. 

3.2.  Related  Works 

Numerous  approaches  for  segmentation  exist  in  the  medical  imaging  literature. 
The  effectiveness  of  these  methods  vary  depending  on  the  application,  noise-level 
and  imaging  modality.  Two  commonly  used  methods  are  deformable  models  and 
clustering-based  methods. 

3.2.1  Clustering  Methods 

Clustering  methods  [5]  are  pattern  recognition  techniques  that  classify  image 
pixels  into  a fixed  set  of  classes  based  on  the  feature-information  present  at  the 
pixel(e.g.,  image  intensity).  Clustering  refers  to  unsupervised  classifier  methods  in 
the  sense  that  no  prior  training  data  is  necessary.  A commonly  used  clustering 
algorithm  is  the  K-means  algorithm  (also  called  ISODATA). 

Given  a metric,  the  K-means  algorithm  iteratively  estimates  a mean  intensity 
for  the  classes  and  allocates  each  pixel  to  the  class  with  mean,  closest  to  its  intensity. 
Even  though  the  K-means  algorithm  is  unsupervised,  it  needs  a good  initialization  for 
the  class-means  to  converge  to  an  optimal  solution.  Further  the  method  is  not  robust 
to  noise  and  intensity  inhomogeneities,  since  it  does  not  impose  any  spatial  or  shape 
constraints  on  the  segmented  regions.  However,  this  lack  of  spatial  modeling  gives  the 
computational  speed  and  ease  of  implementation  usually  associated  with  clustering 
algorithms.  A smoothness  penalty  function  or  a Markov  random  field(MRF)  prior 
([69,  54,  96])  is  often  used  to  enforce  spatial  continuity.  Even  though  these  methods 
can  handle  noise,  intensity  inhomogeneities  and  can  incorporate  prior  information, 
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a major  drawback  is  that  they  do  not  take  into  account  the  shape  of  the  objects 
that  are  to  be  segmented.  State-of-the-art  segmentation  methods  use  clustering  to 
get  a quick  initialization  and  follow  it  with  active  contours  (following  subsection)  to 
capture  the  geometry  of  the  segmentation. 

3.2.2  Active  contour  methods 

Active  contours([52,  87,  42,  10,  44])  are  physically  motivated,  model-based 
techniques  that  can  be  used  to  capture  closed,  lipchitz  boundaries  of  regions  by 
deforming  a curve  under  the  influence  of  internal  and  external  forces.  The  internal 
forces  constrain  the  deforming  contour  to  be  smooth  and  the  external  forces  pull 
the  contour  towards  region  boundaries.  Active  contours  have  the  advantage  that 
they  offer  a coherent  and  consistent  mathematical  description  and  are  robust  to 
noise,  spurious  edges  and  boundary  gaps  due  to  their  incorporation  of  a smoothness 
constraint. 

Active  contour  methods  vary  in  their,  representation  (parametric/implicit)  of 
the  curve  and  modeling(edge-based/region-based)  of  the  external  force.  Parametric 
methods([42,  17,  18])  allow  direct  interaction  with  the  model  and  can  lead  to  a 
compact  representation  for  fast  real-time  implementation.  Adaptation  of  the  model 
topology,  however,  such  as  splitting  or  merging  parts  during  the  deformation,  can  be 
difficult  using  parametric  models.  Implicit  methods([52,  10,  44,  43])  offer  topological 
freedom  for  the  evolving  contour  by  defining  it  implicitly  as  the  0-level  set  of  a higher- 
dimensional function. 

Edge-based,  methods([10,  44,  43,  52])  use  the  image-gradient  information  to 
drive  the  active  contour  towards  image-boundaries.  However,  these  methods  are 
prone  to  “leaks”  in  areas  with  low  contrast,  and  the  solution  is  dependent  on  initial 
contour  placement.  These  limitations  are  overcome  by  Region-based  methods  ([11, 
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84,  97,  80]),  that  use  information  on  regional  statistics  from  the  image.  Many  region- 
based  approaches  are  based  on  the  implementation  of  the  so-called  Mumford-Shah[ M- 
S]  functional [64].  The  problem  is  to  find  a piecewise  smooth  approximation  f for  I 
and  a set  of  boundaries  K that  minimize: 

F(f,  K)  = f (f  — I)2  + a f |Vf|2dx  + £ [ da  (3.1) 

Jn  Jci/k  J K 

where  a and  (3  are  nonnegative  constants  and  JK  dcr  is  the  length  of  K.  In  [11],  Chan 
and  Vese  proposed  an  algorithm  to  treat  a specific  case  of  (3.1),  decomposing  the 
image  into  two  regions  with  common  boundary  C,  and  f constant  in  each  region.  Let 
Qi  and  be  the  regions  inside  and  outside  C.  Let  X\  be  the  characteristic  functions 
for  VL\.  Then  the  following  energy  is  minimized: 

FPC(ci,  c2,  C)  = [ xi(l-ci)2dx+  f X2(i_C2)2dx+  [ da  (3.2) 

Jq  Jci  J c 

where  c,  are  constants.  They  used  a level  set  formulation[65]  to  minimize  the  above 
functional.  Their  use  of  region-based  information  within  an  implicit  framework  gave 
the  advantages  of  segmenting  objects  of  unknown  topology,  robustness  to  noise,  dis- 
continuous edges  and  contour  initialization,  and  detecting  interior  of  objects.  Later, 
in  [87]  they  extend  the  method  to  handle  multiple  regions  and  piecewise  smooth  f. 

However,  the  advantages  of  using  implicit  methods  over  parametric  methods 
come  at  the  cost  of  higher  computational  complexity.  The  embedding  function  for 
the  active  contour  is  usually  constrained  to  be  a signed  distance  function  to  avoid 
numerical  instabilities.  Thus,  usual  implementations  solve  a system  of  two  PDEs. 
One  is  the  equation  for  the  level-set  function  that  advances  the  front.  The  other(re- 
initialization  step)  is  an  initial  value  problem  used  to  preserve  the  signed  distance 
property  of  the  level-set  functions.  Added  to  the  fact  that  there  are  time-step  restric- 
tions on  the  level-set  equation(due  to  the  CFL  condition),  this  turns  out  to  be  time 
consuming  even  in  the  case  of  simple  images  with  distinct  features. 

To  alleviate  these  issues, recently,  variations  to  level-set  implementations  based 
on  the  Mumford-Shah  functional  have  been  proposed([37,  26,  28]). 
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In  [28],  the  image  is  first  preprocessed  using  non-linear  diffusion,  to  remove 
noise  without  losing  image-features.  Then  k-means  is  used  iteratively  with  occasional 
steps  of  mean  curvature  motion,  to  give  a fast  segmentation  of  the  image  into  classes. 

In  [37,  26],  the  characteristic  functions  representing  the  segmentation  are 
found,  using  a phase- fields ([ 63,  3])  based  approach  to  the  Mumford-Shah  functional. 
In  [26],  the  interfaces  are  represented  by  discontinuities  of  piecewise  constant  level 
set  functions  and  each  level  set  function  is  allowed  to  take  only  two  values,  1 or  -1, 
at  convergence.  Their  problem  is  posed  as  the  minimization  of  a smooth  convex 
functional  under  a quadratic  constraint.  Hence  they  avoid  the  re-initialization  step 
used  in  traditional  level  set  implementations.  In  [37],  a phase-field  approximation  for 
(3.2)  is  considered.  Then  using  the  MBO  scheme  of  Osher  et  al. ([61,  60,  62]),  the 
authors  reduce  the  problem  to  simply,  solving  a linear  PDE  with  thresholding. 

In  this  chapter,  we  present  a variational  framework  that  combines  fuzzy  clas- 
sification and  active  contours.  The  fuzzy-membership  functions  and  the  class-means 
form  the  arguments  for  the  energy.  Using  the  membership  functions,  we  define  a hard 
segmentation  of  the  image  which  is  used  to  impose  geometric  constraints  within  the 
energy  functional.  The  membership  functions  and  the  class-means  are  then  computed 
by  numerically  solving  a system  of  PDEs  with  appropriate  initial  and  boundary  con- 
ditions. In  the  case  of  images  with  a high  Signal-to-noise  ratio  (SNR),  some  of  the 
geometric  constraints  can  be  eased  resulting  in  linear  PDE’s  that  are  fast  and  easy 
to  solve.  The  only  initialization  that  is  required  is  for  the  class-means. The  algorithm 
has  shown  a high  robustness  to  noise  and  is  able  to  recover  complex  structures.  We 
have  tested  our  model  on  synthetic  and  MRI  brain  data  and  compare  the  model’s 
performance  with  MRF  based  and  level-sets  based  segmentation  methods. 

Our  work  is  similar  to  [37]  and  [26]  in  our  use  of  a phase-field  like  approach  but 
differs  in  the  sense  that  we  want  a soft  segmentation  of  the  image  into  classes.  Also 
our  method  is  similar  to  region-based  approaches  such  as  ([28,  80])  that  combine 
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fuzzy  clustering  and  level  set  methods.  But  unlike  these  methods,  our  method  is 
variational. 

3.3.  Energy  Functional 

Given  a scalar-valued  image  I : Q — » 3?,  we  look  at  the  problem  of  segmenting 

into  C classes  such  that  I has  a small  intensity- variance  within  each  class.  We 

take  the  soft  segmentation  approach  by  looking  for  C constants  /xk  and  C indicator 

functions  Uk  : — > [0, 1],  (k=l,2,...,C)  such  that: 

c 

£ u*  = 1 <3-3) 
k=l 

Denote  /x  = (/x i,  /x2,  •••,  /x c)  and  U = (ux,  u2, ....,  uc). 

We  look  at  the  following  energy  functional: 

c c 

Ec[/b  U]  = f ^2  M1  - A*k)2  + p'52  I Vuk|2dx  (3.4) 

k=l  k=l 

For  each  x E and  (k— 1,2,..,C),  the  first  term  pushes  Uk(x),  towards  0 if  I(x)  differs 
from  /Xk  and  towards  l(due  to  the  constraint  (3.3))  if  I(x)  is  close  to  /xk.  The  second 
term  in  the  energy  acts  as  a regularizing  term  for  {uk}k=1,  P acts  as  a balancing 
parameter  between  the  two  terms.  The  constraint  (3.3)  prevents  the  minimizer  from 
taking  a trivial  solution,  i.e  U = 0. 

In  order  to  decrease  the  number  of  arguments  in  (3.4)  and  also  to  remove  the 
constraint  (3.3)),  we  instead  assume  C to  be  a power  of  2 and  look  for  N = log2C 
class-functions  Lk  : f2  — > [0,  |]  that  minimize  the  following  energy: 

» C N 

Ec [/x.,  L]  = / ^Ak(I-/Xk)2+/3^|VLk|2dx  (3.5) 

*'n  i=l  k=l 

where  L = (Li,  L2,  ...LN)  and  Ak  are  defined  as  follows. 

If  we  let  B;  be  the  N-bit  binary  representation  for  (i-1),  then  we  can  define 
Bik  to  be  the  kth  bit  of  B;  and  A,  = JIk=i 

where  Mjk  = cos(Lk) (respectively,  sin(Lk))  if  Bik  =0 (respectively,  1). 

It  can  be  seen  that  {A;}^  act  as  the  indicator  functions  and  satisfy  XliLi  = 
1.  To  illustrate,  in  the  2-phase  case,  the  energy  E2  is  given  by: 
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E2[a*i,  Ai2,  Li]  = / cos(L1)2(I-/ii)2  + sin(L1)2(I-/x2)2  + /5|VL1|2dx 
Jq 

In  the  4-phase  case,  the  energy  E4  is  given  by: 

E4|/b,  A*2,  ^3, A4, Ll5  L2]  = / Ai(I  — fix)2  + A2(I  — /i2)2  + A3(I  — //3)2  -I-  A4(I  — pf)2  + 

Jn 

/0(|  VL4 12  + |VL2|2)dx 

where, 

A4  = cos(Li)2cos(L2)2 
A2  = cos(Li)2sin(L2)2 
A3  = sin(L!)2cos(L2)2 
A4  = sin(Li)2sin(L2)2 

3.4.  Phase  Functions  and  Geometric  Constraints 
From  the  indicator  functions  {Ak}^=1  in  (3.5),  we  define  the  N phase-functions 
as  follows: 

$1  = cos(Lj)2  — sin(Lj)2 

for  (i=l,2,...N) 

Using  these  functions,  we  can  apply  geometric  constraints  on  the  segmentation  {fhJkLi 
such  as  area,  perimeter  and  curvature,  within  the  energy  (3.5).  We  demonstrate  this 
in  the  4-phase  case: 

To  start  with,  we  define  our  hard  segmentation  {r2k}^_1  using  the  maximum 
membership  value: 

fik  = {x|Ak  > A;|k^iil<i<4} 

To  simplify  computations  we  use  the  phase  functions  <f>i  and  $2  to  represent  {Ok}k=i 
(Fig.3.1).  It  is  straightforward  to  show  that: 
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fii  = {x|A4  > Ai  and  A4  > A2  and  A4  > A3}  = {x|$i  < 0 and  $2  < 0} 

Fl2  = {x|A2  > Ai  and  A2  > A3  and  A2  > A4}  = {x|$i  > 0 and  $2  < 0}  (3.6) 

— {x|A3  > Ai  and  A3  > A2  and  A3  > A4}  = {x|$4  < 0 and  <f>2  > 0}  (3.7) 

Q4  = {x|A4  > A2  and  A4  > A3  and  Ai  > A4}  = {x|$4  > 0 and  $2  > 0}  (3.8) 

The  functions  $;  are  the  soft-analogies  to  the  level-set  functions  used  in  Chan-Vese’s 
work  on  multiphase  segmentation  [87]. 

For  (k=  1 ,2,3,4) , let  Sf2k  denote  the  boundary  of  Fl^.  Its  perimeter  (lk)  and 
area(Sk)  can  be  easily  computed.  For  instance,  for  k=l: 

U = [ |V(H($1)H($2))|dx 
Jn 

Si  = [ H($1)H(«F2)dx  (3.9) 

in 

where  H(t)  is  the  heaviside  function. 

Let  T = Uk=i  denote  the  edge-set  between  the  classes. 

We  can  use  the  following  term  to  constrain  the  length  of  F : 

lr=  [ |VH($1)|  + |VH($2)|dx 
in 

It  is  to  be  noted  that  the  above  expression  does  not  calculate  the  actual  length  of  F 
since  it  weighs  some  parts  of  T more  than  others,  as  mentioned  in  [37].  An  alternative 
way  can  be  formulated  as  in  [37]: 


k=l 


In  a similar  way,  the  above  computations  can  be  extended  for  the  C-phase  case. 
Further,  if  one  wants  completely  stay  in  the  fuzzy  setting,  constraints  from  fuzzy 
geometry[71]  such  as  area,  perimeter  can  be  used. 

3.5.  Euler-Lagrange  Equations 

For  the  sake  of  simplicity  of  presentation,  we  demonstrate  the  minimization 
problem  for  the  4-phase  case  only.  The  functional  (3.5)  with  a length  constraint  on 
the  edge  set  F is: 
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Figure  3.1:  Phases  Qk  and  edge  set  F (union  of  all 

curves  seen  in  the  figure) 


E4[^i,  M2)  1^3,  A4,  Li,  L2]  — / Ai(I  — //i)2  + A2(I  — /i2)2  + A3(I  — //3)2  + A4(I  — /^4)2  + 

J n 

/5(|VL4|2  + |VL2|2)  + A(|VH($X)|  + |VH($2)|)dx 
where  A is  a balancing  parameter  and  Aj  are  as  defined  in  (3.5).  The  corresponding 
Euler-Lagrange  equations  (necessary  condition  of  optimality)  for  E4[/i,  L]  are 
/?ALi  = sin(2L!)[cos2(L2)(I  - /i3)2  + sin2(L2)(I  - /x4)2  - cos2(L2)(I  - fii)2 

- sin2(L2)(I  - /x2)2]  - A5(4>1)div(^j|i]j),  (3.10) 

/3AL2  = sin(2L2)[cos2(Li)(I  - //2)2  + sin2  (Li)  (I  - /x4)2  - cos2(Li)(I  - /xi)2(3.11) 

- sin2(L!)(I  - /x3)2]  - A<S($2)div(||[!^),  (3.12) 

(3.13) 


dU 

dn 


dh2 

dn 


= 0 on  dtt, 


with 
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and 


/n  Ai  • I dx  _ fn  A2  • I dx  fQA 3-I  dx 


u — --  ..  _ ..  _ Jq  -1 

A*1  — r a > > ^2  — — f — x — 7~  > 4*3  — — r — a — T- > A4  — 


In  Ai  dx 


In  A2  dx 


In  A3  dx 


In  A4 


where  5 is  the  dirac  delta  function.  In  our  implementation,  we  use  a smooth  version 
of  £(*),  given  by  <$e(t)  = -(1glVt)i- 

3.6.  Numerical  Scheme 

Note  that  the  Euler-Lagrange  equations  for  Li  and  L2  are  of  the  type 

r\ 

f(u)  - Au  = 0,  — ' - = 0 on  dQ. 

on 

The  corresponding  gradient  flow  is  given  by 

ut(t)  = 


f(u)(t)  + Au(t),  = 0. 


This  equation  is  integrated  using  a semi-implicit  scheme.  Let  uy  denote  the  value  of 
u at  the  pixel  xs  = ih,  yj  = jh,  where  h is  the  pixel  size.  Denote  u(xj,yj,tn)  by  u^. 
We  use  the  following  scheme  to  approximate  ufd1: 

<? 1 = - MMj)  - + u^+1  + - 4<+‘)), 

<t‘(l  + 4A t)  = ufj  - A «(/«,-)  - «+1J  + uSLy  4-  «?J+1  + (3.14) 

where  At  is  the  time-step. 


Once  every  few  iterations(typically  50),  the  class-functions  {Lfc}f=1  need  to 
be  re-scaled  to  let  each  of  the  indicator  functions  {AJ-Ij,  take  the  range  [0,1].  This 
rescaling  prevents  the  functions  from  becoming  flat  especially  in  the  case  of  noisy 
and  low-contrast  images. 

This  is  similar  to  the  re-initialization  step  in  many  level-set  methods[87,  11] 
that  require  the  level-set  functions  to  be  signed  distance  functions  close  to  the  evolving 
front.  It  is  to  be  noted  that,  as  of  now  we  have  no  theoretical  justification  for  this 
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Figure  3.2:  I)Noisy  image  II)  The  white  curve  shows 

the  edges. 


re-scaling  step,  but  have  found  it  to  be  necessary  when  handling  noisy  images  as  in 
Fig.3.4(IV). 

We  also  wish  to  comment  that,  during  minimization,  the  active  contour  cannot 
be  explicitly  determined  from  the  phase  functions  <Fj,  since  |V<Fj|  are  not  constrained. 
Hence,  the  geometric  constraints  introduced  within  the  functional,  influence,  not  only 
the  active  contour,  but  also  a neighborhood  of  level  curves  around  the  contour.  The 
width  of  this  neighborhood  depends  on  the  level  of  uncertainty  that  exists  for  pixels 
in  the  image.  A good  indicator  for  the  uncertainty  in  the  object  boundaries  is  the 
width  of  the  support  of  the  functions  6e(4>j) (defined  in  the  previous  section). 

3.7.  Experimental  Results 

This  section  shows  the  results  of  applying  our  fuzzy  segmentation  model  to 
synthetic  and  MRI  data.  The  synthetic  experiments  demonstrate  the  model’s  capa- 
bility to  constrain  the  geometry  of  the  segmented  objects.  Hence,  similar  to  active 
contour  methods,  the  model  has  shown  a high  tolerance  to  noise(Fig.3.2,  Fig. 3. 4) 
and  can  be  used  group  objects(Fig.3.4).  We  also  validate  the  model’s  performance 
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Figure  3.3:  I)  European  lights  image  II)  Seg- 
mentation obtained  with  A = 10,  /3  = 3. 


on  simulated  MR  images  taken  from  the  Brainweb  database  [20],  for  varying  levels  of 
noise  and  intensity  inhomogenietes.  Finally,  we  present  preliminary  results  on  a real, 
3D  MRI  dataset. (Fig. 3. 9) 

The  number  of  classes  C is  assumed  to  be  known.  Like  most  clustering  meth- 
ods, the  algorithm  needs  an  initial  guess  for  the  class-means  and  within  our  ex- 
periments, the  results  have  been  found  to  be  fairly  robust  to  the  initial  guess.  In  our 
implementation,  we  compute  the  histogram  of  the  given  image  I with  C bins  and  set 
Hk  equal  to  the  center  of  the  kth  bin.  The  class-functions  {L*,}^  are  usually  picked 
to  be  | for  all  images  and  the  parameters  /3  and  A need  to  be  tuned  for  each  image. 
These  parameters  have  to  be  carefully  picked  especially  for  noisy  images,  since  a very 
large  value  can  result  in  loss  of  edge  features  and  likewise  a small  value  can  lead  to 
a noisy  segmentation.  In  this  work,  we  estimate  j3  and  A on  a trial  and  error  basis. 

Also,  convergence  can  be  considered  to  be  obtained  when  the  change  in  the 
objective  function  is  less  than  some  threshold,  or  when  the  change  in  the  indicator- 
function  values  is  less  than  some  threshold.  The  latter  criterion  is  used  in  this  chapter. 
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Figure  3.4:  I)Original  image.  II)  A = 0, 

P = 1 HI)  A = l,p  = 3 IV)  A = 20./3  = 1 

3.7.1  Synthetic  Images 

In  Fig.3.2,  we  want  to  segment  a noisy  image(with  additive  gaussian  noise[0 
mean, .4  standard  deviation])  into  two  classes.  Starting  with  initial  values  for  the  two 
class  means,  computed  as  indicated  above,  we  obtain  a segmentation(edges  shown  by 
white  curve  in  II).  For  this  image  we  chose  A = 5,  P — .5  as  the  parameter  values  and 
a time  step  At  = 2.  We  obtained  stable  convergence  in  300  iterations.  It  is  to  be 
noted  that  without  the  length  term  (i.e.  when  A = 0),  p has  to  be  weighted  more  to 
handle  the  noise  and  some  of  the  sharp  corners  cannot  be  captured  in  the  resulting 
segmentation.  Further,  use  of  the  length  term  allows  the  segmented  objects  to  have 
smooth  boundaries. 
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Figure  3.5:  I)Original  image.  II)  A = 0, 

0 = 1 III)  A = 1,0  = 3 IV)  A = 20,0  = 1 


Fig. 3. 3 and  Fig. 3. 4 demonstrate  the  ability  of  the  model  to  group  objects  like 
active-contour  methods.  In  Fig. 3. 3,  we  see  a segmentation  for  the  ’european  lights’ 
image  with  A = 10  and  0 = 3.  Use  of  the  length  term  resulted  in  large  connected 
regions  being  preferred  in  the  segmentation. 

In  Fig.3.4.  (II  and  III)  ,we  obtain  different  segmentations  by  varying  the  para- 
meters A and  0.  In  (III),  we  see  that  for  the  value  of  A and  0 chosen,  the  segmentation 
has  captured  gaps  in  the  boundary  however  at  the  cost  of  losing  some  corners.  In 
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Figure  3.6:  Ground  truth  for  fuzzy  classes  I)  CSF  II)  GM  III)WM 

(IV),  the  function  5e(<I>)  is  displayed  as  a grayscale  image.  The  width  of  the  white 
band  indicates  the  fuzziness  of  the  object  boundary. 

The  importance  of  the  re-scaling  step  can  be  seen  in  Fig.3.5.  In  the  case 
of  an  image  where  the  class-means  are  very  close  to  each  other (e.g.,  with  noisy  or 
low-contrast  images,  the  phase-function  tends  to  flatten  out  due  to  the  smoothing 
terms.  The  re-scaling  step  is  used  to  counter  this  problem.  In  Fig. 3. 5(1),  a very 
noisy  version  (speckled  noise  75%)  of  Fig.3.4(I)  is  shown.  With  the  rescaling  step,  a 
reasonable  segmentation  is  obtained  in  Fig.3.5(II(a)).  The  fuzziness  of  the  segmented 
boundary  is  seen  in  11(b).  In  contrast,  without  the  rescaling  step,  the  phase  function 
has  flattened  out(III(b))  and  the  segmentation(III(a))  computed  is  not  correct. 

3.7.2  MRI  Brain  Data 

The  study  of  many  brain  disorders  involves  accurate  tissue  segmentation  from 
magnetic  resonance  (MR)  images  of  the  brain.  Automated  and  reliable  segmentation 
of  MR  images  is  complicated  due  to  partial  volume  effects,  noise  and  a spatially 
varying  intensity  inhomogeneity.  In  this  section,  we  validate  the  application  of  our 
model  to  automatic  MRI  segmentation. 

Since  the  ground  truth  of  segmentation  for  real  MR  images  is  not  usually 
available,  it  is  impossible  to  evaluate  the  segmentation  performance  quantitatively, 
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but  only  visually.  However,  Brainweb  provides  a simulated  brain  database  (SBD) 
including  a set  of  realistic  MRI  data  volumes  produced  by  an  MRI  simulator  from 
known  (ground  truth)  models.  We  evaluate  the  performance  of  our  segmentation  on 
some  brainweb  data  by  comparing  it  with  the  given  ground  truth.  We  also  present 
comparisons  with  some  current  fuzzy  segmentation  methods. 

We  picked  the  64th  transversal  slice  from  a simulated  3D  volume  with  di- 
mensions (181  x 217  x 181)  and  segmented  it  into  4 tissue  classes(background(BG), 
cerebrospinal  fluid(CSF),  grey  matter(GM)  and  white  matter(WM)). 

Fig. 3. 7 shows  the  segmentation  that  was  obtained.  (I)  and  (II)  are  the  results 
for  the  image  with  no  inhomogeneities,  for  3%  and  9%  level  gaussian  noise,  and  (III) 
is  the  result  for  9%  noise  and  40%  intensity  inhomogeniety.  The  results  in  (I)  and 
(II)  seem  accurate,  atleast  visually.  We  will  validate  these  results  in  the  following 
discussion.  In  (III),  some  pixels  in  the  lower  part  of  the  brain  have  been  misclassified, 
due  to  large  intensity  inhomogenities  in  that  region.  It  is  to  be  noted  that  our  model 
cannot  handle  large  inhomogenieties.  However,  it  is  possible  to  extend  our  model  to 
handle  inhomogenieties,  by  modifying  the  objective  functional(3.4),  similar  to  Pham 
et.al([70]).  Neglecting  the  BG  class,  we  compare  the  fuzzy  membership  functions(for 
CSF,  GM,  WM)  computed  with  our  model,  with  the  available  ground  truth  for  the 
fuzzy  classes.  (Fig. 3. 6). 

The  fuzzy  segmentation  for  a class  k is  measured  by  the  quantity  £k  defined 
by: 

f = Exes  K(x)  -gk(x)[ 

Card(S) 

where  S is  the  set  of  brain  pixels  (we  get  this  set  after  applying  a brain  mask  to  the 
image)  , Uk  is  the  membership  function  computed  by  the  model  for  the  kth  class  and 
gk  is  the  corresponding  ground  truth.  Here,  the  classes  k=l,2,3  correspond  to  the 
tissues  CSF,  GM  and  WM  respectively. 
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Table  3.1:  Comparisons  of  our  method  with  FCM  and  FMS  for  3%, 5%, 7%  noise  levels 


Tissue/Noise 

3% 

FCM 

5% 

7% 

3% 

FMS 

5% 

7% 

3% 

FAC 

5% 

7% 

CSF(&%) 

2.72 

3.6 

4.50 

2.18 

2.56 

3.89 

8.29 

8.69 

9.05 

GM  (fc%) 

7.54 

10.97 

15.14 

7.23 

9.28 

12.41 

5.72 

5.78 

5.98 

WM(&%) 

5.33 

8.07 

11.22 

4.70 

5.97 

9.30 

3.76 

3.96 

4.30 

Table  3.2:  Effect  on  average  error,  £,  for  3%,  5%,  7%,  9%  noise  levels 


Length 
term,  A 

3% 

5% 

7% 

9% 

0 

5.92 

6.14 

6.44 

6.88 

.2 

6.02 

6.20 

6.47 

6.90 

We  computed  the  fuzzy  classification  error  for  different  noise  levels  and  com- 
pared the  error  obtained  by  FCM  and  Fuzzy  Markovian  segmentation  method(FMS)[73] 
(Table.3.1).  Our  model  outperformed  FCM  and  FMS  in  its  ability  to  handle  noise, 
which  is  seen  from  the  small  changes  in  that  occurs  with  increasing  noise  levels. 
Further,  the  errors  for  WM  and  GM  are  much  lesser  in  our  case.  On  the  downside, 
our  smoothing  constraints  affects  the  recovery  of  small  structures  such  as  in  the  CSF 
class,  which  is  evident  from  our  error^). 

In  Fig.3.8  for  the  same  image  with  9%  noise,  addition  of  a length  constraint  has 
given  a segmentation  with  smoother  boundaries.  Table.3.2  shows  the  effect  of  length 
term  on  the  average  error,  f = mean{£ i,f2,£3}-  The  table  shows  that  the  addition 
of  the  length  term,  causes  small  distortions  in  the  fuzzy  membership  functions.  But 
these  distortions  become  negligible  for  large  noise  levels.  Moreover,  the  average  error 
of  our  model  with  the  length  term,  is  lesser  than  that  of  FCM  and  FMS.  Hence  it  is 
viable  to  use  the  length  term  to  extract  smooth  boundaries. 

In  Fig.3.9,  a 3D  visualization  is  shown  of  the  cortex  which  corresponds  to 
one  of  the  phases  computed  from  a 4-phase  segmentation  on  a 3D-T1  MRI  data. 
The  3D  data  was  obtained  courtesy  Siemens  Corporate  Research,  New  Jersey.  In 
this  experiment,  the  smoothing  parameter  j3  was  given  a large  weight,  to  prevent 
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pixels  belonging  to  the  skull(that  have  an  intensity  profile  similar  to  WM),  being 
misclassified  into  the  WM  phase.  As  a result,  the  rendered  surface  has  not  captured 
many  of  the  cortex’s  features.  A common  way  to  handle  this  problem  is  to  strip 
the  skull-pixels  from  the  data[78],  and  then  apply  the  segmentation  algorithm.  We 
expect  a better  3D  MRI  brain  segmentation  from  our  model  after  the  skull  stripping 
step. 

3.8.  Conclusion 

We  looked  at  a variational  method  for  unsupervised  classification  of  a given 
image  into  a fixed  number  of  classes.  The  resulting  segmentation  is  fuzzy,  in  the  sense 
that  it  is  characterized  by  membership  functions  for  each  class.  These  membership 
functions  give  a hard  segmentation  that  is  used  to  impose  geometric  constraints  on 
the  shape  of  the  classes  within  the  energy  framework.  These  constraints  give  the 
robustness  to  noise  that  is  usually  associated  with  active-contour  methods.  Use  of 
a clustering  based  framework  allows  use  of  region-based  information,  resulting  in 
fast  and  stable  convergence  to  a global  minimum.  State-of-the-art  active  contour 
methods  use  FCM  as  a pre-processing  step  to  get  a closer  initialization  to  the  global 
minimum.  Our  model  combines  FCM  and  active  contours  in  a unified  framework  with 
a simple  implementation  and  gives  results  comparable  to  traditional  clustering-based 
and  active-contour  methods. 

Our  current  work  involves  adapting  this  model  to  extract  local  structures  and 
will  form  part  of  future  communication  on  this  work.  We  are  also  interested  to  look 
at  the  existence  and  uniqueness  properties  of  the  minimizers  of  energy  (3.5)  and 
investigate  the  significance  of  the  re-scaling  step  of  the  class-functions. 
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Figure  3.7:  I)Original  image  with  40%  intensity  inhomogeniety  II)Segmentation 

III)Original  image  with  40%  intensity  inhomogeniety  and  gaussian  noise. 
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Figure  3.8:  I)Image  with  9%  noise  level  II)Segmentation  with  no  length  term(rough 
boundaries)  III) A = .2(smooth  boundaries) 


Figure  3.9:  3D  visualization  of  the  boundary(  cortex)  of 
the  4th  phase  after  a 4-Phase  segmentation. 


CHAPTER  4 

JOINT  SEGMENTATION  AND  REGISTRATION  OF  FMR  IMAGES 

Motion  can  significantly  impact  signal  changes  in  functional  MR  images(fMRI), 
and  affect  the  detection  of  task  activation.  To  minimize  motion-induced  noise,  we 
propose  a feature-based  image  registration  model.  Due  to  T2  weighted  signal  loss,  de- 
creased resolution,  and  low  contrast  in  fMR  images,  the  feature  used  for  registration 
cannot  be  detected  reliably  by  using  only  image  information.  In  many  medical- 
imaging applications,  information  on  the  shape  of  the  feature  to  be  detected  is  avail- 
able. We  use  this  information  with  image-intensity  to  reliably  segment  the  required 
feature,  in  a low-resolution  image.  The  model  is  formulated  as  the  minimization  of 
an  energy  functional  that  depends  on  the  image  gradient  and  the  given  shape,  so 
that  the  boundary  of  the  object  captured  occurs  at  high  gradients  and  is  as  close  as 
possible  to  the  given  shape.  Further,  the  transformation  that  aligns  the  segmented 
feature  with  the  prior-shape  is  computed  as  part  of  the  model  and  is  used  in  the 
alignment  of  a time  series  of  images.  The  model  was  tested  on  both,  synthetic  data 
and  fMRI  brain  data,  with  good  results. 

4.1.  Introduction 

Functional  MRI  is  an  extension  of  MRI  that  measures  changes  in  cerebral 
blood  flow  and  oxygen  level  that  reflect  localized  changes  in  brain  activity  induced 
by  sensory,  motor,  or  cognitive  tasks,  by  acquiring  a rapid  series  of  MRI  images  of 
a subject’s  brain  under  different  task  conditions.  Due  to  the  long  scanning  time 
for  completing  the  series  of  images,  head  movement  is  unavoidable.  One  of  the 
challenges  facing  fMR  imaging  is  to  distinguish  between  task-activation  induced  and 
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head-motion  induced  signal  changes.  Since  task-induced  signal  changes  by  only  5- 
10%  above  the  mean  signal  intensity  in  fMRI,  even  small  motions  can  severely  impede 
optimal  detection. 

To  minimize  the  effects  of  motion-induced  noise,  the  series  of  fMR  images  have 
to  be  first,  registered,/ aligned.  There  have  been  a number  of  registration  methods 
developed  for  medical  images.  For  recent  surveys,  see  ([35,  51]).  Generally  speaking, 
these  methods  are  classified  into  intensity  based  and  feature  based  methods. 

Intensity  based  registration  methods  ([39,  92],  and  references  therein)  align  a 
pair  of  images  by  maximizing  intensity  correlation  or  by  minimizing  intensity  vari- 
ance. These  methods  rely  on  the  implicit  assumption  that  when  two  images  are 
perfectly  aligned,  their  voxel  by  voxel  difference  is  very  small.  It  has  been  pointed 
out  in  [6]  that  this  assumption  is  not  strictly  valid  for  fMR  images  in  either  task 
activation  or  in  resting  state  when  pairwise  image  comparisons  are  made  within  the 
same  time  series  of  images  or  across  different  functional  image  sets. 

A more  general  and  popular  method  for  intensity  based  registration,  is  maxi- 
mizing mutual  information  of  the  reference  and  the  transformed  new  image  ([92,  22, 
34,  82,  88,  90,  91]).  However,  due  to  the  complex  nature  of  intensity  variations (e.g., 
differences  in  signal  intensity-distributions  during  task  activation  and  a rest  scan,  and 
intensity-inhomogeneities),  maximizing  the  mutual  information  between  two  fMR  im- 
ages does  not  give  a good  registration. 

Feature  based  methods  first  segment  important  features  in  the  images  and 
then  match  these  features  in  different  images  by  adjusting  the  motion  parameters 
([6,  47,  85]  and  references  therein).  Feature-based  methods  are  not  adversely  affected 
by  intensity  complexity  and  can  be  used  for  registering  images  of  different  modality. 
However,  the  quality  of  registration  depends  on  the  segmentation  step,  and  for  im- 
ages such  as  fMRI,  feature-segmentation  is  challenging.  Image  segmentation  can  be 
broadly  classified  into  edge-based  and  region-based  methods. 
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Edge-based  segmentation  methods  use  the  information  of  sharp  intensity  changes. 
([10,  44,  43,  52]).  Region-based  methods  make  use  of  homogeneity  on  the  statistics 
of  the  regions  being  segmented  ([11,  64,  84,  97]).  Models  that  integrate  edge  and 
region-based  methods  have  also  been  developed.  ([66,  68]). 

But  for  fMRI,  due  to  T2  weighted  signal  loss  and  decreased  resolution,  the 
boundaries  of  the  features  are  both  low  contrast  and  low- resolution.  Since  edge-based 
models  often  ‘leak’  through  gaps(points  of  low  image  gradients)  in  the  boundary, 
and  region-based  methods  rely  on  the  homogeneity  of  the  regions  separated  by  the 
boundary,  both  methods  cannot  detect  features  in  fMR  images  reliably. 

In  some  medical  image  processing  problems,  the  information  about  the  shape 
that  we  want  to  detect  is  available,  which  motivates  us  to  look  for  a method  that 
incorporates  prior  shape  information  into  the  process  of  segmentation. 

Several  methods  for  incorporating  prior  shape  information  into  boundary  de- 
termination have  been  developed  (e.g.,  see  [21,  46,  79,  89,  95]).  In  [21]  and  [89],  a 
statistical  model  of  shape  variation  is  constructed  from  a set  of  corresponding  points 
across  the  training  images,  and  then,  a Bayesian  formulation  based  on  this  prior 
knowledge  and  the  edge  information  of  the  image  is  employed  to  find  the  object 
boundary.  In  [79]  an  elliptic  Fourier  decomposition  of  the  boundary  is  used  to  in- 
corporate the  global  shape  information  into  the  segmentation.  In  [95],  a method  of 
deformable  templates  is  proposed  for  feature  extraction  from  faces.  In  this  model,  a 
parameterized  template  describes  the  feature  of  the  interest.  The  deformable  tem- 
plate interacts  dynamically  with  the  image  by  varying  its  parameter  values  to  mini- 
mize the  energy  function,  which  links  edges,  peaks,  valleys  in  the  image  intensity  to 
the  corresponding  properties  of  the  template. 

Recently,  statistical  shape  knowledge  has  been  incorporated  into  edge  based 
or  region  based  active  contours.  In  [46],  Leventon,  Grimson  and  Faugeras  extended 
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geometric  active  contours  ([10,  44])  by  incorporating  shape  information  into  the  evo- 
lution process.  They  built  a statistical  shape  model  from  a training  set  of  curves 
represented  by  signed  distance  functions.  The  segmentation  process  embeds  an  ini- 
tial curve  as  a level  set  of  a higher  dimensional  surface,  and  evolves  the  surface  based 
on  image  gradients  and  curvature,  and  towards  a maximum  a posteriori  estimate  of 
shape  and  pose.  In  [22]  Cremers,  Schnrr  and  Weickert  incorporated  statistical  shape 
knowledge  into  the  Mumford-Shah  segmentation  scheme  [64]  by  minimizing  a func- 
tional that  includes  the  Mumford-Shah  energy  and  the  shape  energy  corresponding 
to  the  Gaussian  probability. 

Prior  statistical  shape  information  is  also  used  to  control  a total  deformation 
of  the  active  contour  formed  by  global  transformation  (e.g.,  rotation,  scaling,  and 
translation)  plus  a local  deformation  (see  [72,  75]). 

Our  basic  idea  is  also  to  influence  a geometric  active  contour  with  a vector  field 
that  depends  on  shape  priors.  Our  model  differs  in  that,  we  use  a variational  approach 
instead  of  a probability  approach.  This  makes  it  possible  to  discuss  the  existence 
of  the  model  solution  in  BV  (bounded  variation)  space.  This  chapter  presents  a 
new  variational  framework  that  combines  shape  information  with  an  edge-/region- 
based  term,  offering  a meaningful  way  to  capture  the  boundaries  for  low-contrast /- 
resolution  images.  Given  an  image  I and  prior-shape  S,  the  problem  can  be  stated  as 
finding  an  edge  C in  I and  a transformation  T defined  on  C,  such  that  the  transformed 
edge  ToC  is  close  to  S.  Typically,  formulated  as  an  energy  functional  minimized  over 
C and  T,  such  a model  would  have  two  terms: 

A)  An  edge  term  P that  captures  object  boundaries  regions  in  a given  image 

B)  A shape-similarity  term  Dg  to  measure  the  closeness  of  the  segmented  feature  to 
the  given  shape. 

This  functional  can  be  written  as 


E[C,T]  =P(C;I)  + A Dg(ToC) 
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where  A balances  the  edge  and  shape-similarity  terms.  In  our  work: 

A)  The  shape-similarity  term  Ds  = d§,  where  ds  is  the  distance  function  for  S. 

B)  The  transformation  T is  rigid. 

C)  We  consider  two  models  for  the  edge-term  P,  one  is  edge-based  (Geometric  active 
contours^ GAC),[44])  and  the  other,  region-based(Mumford-Shah[M-S ] functional, [64]). 

The  segmented  curve  C is  driven  towards  boundaries  of  the  image  I with  its 
shape  being  constrained  to  be  close  to  the  prior  shape  S,  which  is  particularly  useful 
in  completion  of  missing  image-boundaries(such  as  in  fMR  images). 

In  the  following  sections,  we  present  our  two  models.  One  uses  the  edge-based 
GAC  and  the  other  uses  the  region-based  M-S  functional.  We  discuss  the  formulation, 
implementation  and  experimental  results  for  each  of  the  models. 


Let  I : Q 5i,  Q C K be  an  image.  Let  I*  is  a gaussian  smoothed  version 
of  the  image  I.  Let  C(p)  = (x(p),  y(p))  (p  G [0,1])  be  a differentiable  parameterized 
curve  in  fh  The  geometric  active  contour  model  is  formulated  as  the  minimization 
of  the  energy  function: 


where,  the  integral  is  with  respect  to  the  arc  length  function  of  C and 


The  above  functional  is  minimized  when  |VI*|  takes  large  values  on  C.  The  main 
advantages  of  GAC  is  its  ability  to  directly  generate  closed  parametric  curves  from 
images  and  the  incorporation  of  a smoothness  constraint  that  provides  robustness  to 
noise  and  spurious  edges.  The  main  disadvantage  is  its  dependence  on  the  image- 
gradient  which  leads  to  ‘leaks’  across  image  boundaries  in  the  case  of  low-contrast /- 
resolution  images(such  as  an  fMR  image). 


4.2.  GAC  based  segmentation  with  shape 
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Let  S be  a curve,  representing  the  shape  we  expect  of  the  boundary.  In  order 
to  add  the  shape  term  to  the  GAC,  we  introduce  the  energy  functional  E(C,  /x,  R,  T), 
where  (/x,  R,  T)  is  a rigid  transformation  of  C through  scaling  /x,  rotation  R and 
translation  T: 

E(C,  R,  T)  = £ g(|  Vr(C)|)  + ^d2(^RC  + T)  ds  (4.1) 

where,  A > 0 is  a parameter,  and  d(x,  y)  is  the  distance  of  the  point  (x,y)  from  S.  In 
the  numerical  computation,  this  distance  function  is  obtained  by  the  fast  marching 
method  proposed  by  Sethian.  The  first  term  on  the  right  hand  side  of  (4.1)  comes 
from  GAC  and  measures  the  amount  of  high  gradient  under  the  trace  of  the  curve 
C.  The  second  term  measures  the  closeness  of  C to  the  prior  shape  S. 

The  curve  C and  the  transformation  parameters  (/x,  R,  T)  evolve  to  minimize 
E(C,/x,R,  T).  The  evolution  is  done  by  gradient  descent  with  respect  to  (/x,  R,  T) 
and  first  variation  descent  with  respect  to  C. 

The  gradient/variation  descent  for  the  energy  functional  (4.1)  is  performed  by 
the  following  system: 


Ct  = -v. n,  C(0,p)  = Co(p), 

(4.2) 

/Xt  = — \ j)  dVd  • RC  ds,  n( 0)  = /x0, 

(4.3) 

>t  = — A/x^dVd-^C  ds,  0(0)  = 0O, 

(4.4) 

Tt  = —A  dVd  ds,  T(0)  = T0, 

(4.5) 

where  n is  the  outward  unit  normal  to  C, 


v = Vg  • n T gk  + A/x(d  • Vd)  • (Rn)  + Ad2k, 
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and  k is  the  curvature  of  the  curve  C.  In  Eqs.(4.3)  - (4.5), 

the  function  d is  evaluated  at  //RC(p)  + T , while  the  functions  g,  k,  and  n are 
evaluated  at  C(p). 

We  then  use  a finite  difference  scheme  to  numerically  approximate  the  above 
system.  In  problems  of  curve  evolution  the  level  set  methods  have  been  used  exten- 
sively, because  it  allows  for  cusp,  corners,  and  automatic  topological  changes.  In  the 
following  discussion,  we  consider  the  level  set  formulation  for  the  model  (4.1)  and 
give  a proof  for  the  existence  of  minimizers. 

4.2.1  Level  Set  Formulation 

Represent  a contour  C by  the  zero  level  set  of  a Lipschitz  function  u such 
that  (x|u(x)  > 0}  is  the  set  inside  C.  Let  H(z)  be  the  Heaviside  function,  that  is 
H(z)=l  if  z > 0,  and  H(z)=0  if  z < 0,  and  6 = H'(z)  (in  the  sense  of  distribution) 
be  the  Dirac  measure.  Then,  the  length  of  the  zero  level  set  of  u in  the  confor- 
mal metric  ds  = g|C'(p)|dp  can  be  computed  by  /ng|VH(u)|  = Jn5(u)g|Du|.  The 
similarity  of  the  shapes  between  the  zero  level  set  of  u and  C*  can  be  evaluated  by 
fQ  S( u)d2(//Rx  + T))dx,  where  the  function  d is  the  same  as  in  (4.1).  Therefore,  the 
level  set  formulation  of  our  variational  approach  can  be  formulated  by: 


(4.6) 


The  evolution  equations  related  to  the  Euler-Lagrange  equations  for  (4.6)  are 


— = 0,  x 6 dfl,  t > 0, 


u(x,  0)  = u0(x),  x E Q, 


(4.7) 
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0/1  /»  1 T~) 

^ = -Ay  5(u)//dVd  • (— x)|Vu|dx,  t>0,  9(0)  = 90,  (4.9) 

dT  r 

— = -A  / <5(u)dVd|Vu|dx,  t > 0,  T(0)  = T0.  (4.10) 

d’t  J Q 

where  R is  the  rotation  matrix  in  terms  of  the  angle  9,  and  the  function  d is  evaluated 
at  //Rx  4-  T. 

We  use  a finite  difference  scheme  to  numerically  approximate  the  above  sys- 
tem, detailed  in  our  paper([12]). 

4.2.2  Proof  of  Existence 

Here,  we  discuss  the  existence  of  the  solution  to  the  minimization  in  equation 
(4.6).  The  same  argument  can  also  be  applied  to  equation  (4.1).  As  suggested  in  [11] 
we  may  rewrite  (4.6)  in  the  form 

min  f {g(|  VI| (x)  + ^d2(/iRx  + T)}|DyE|,  (4.11) 

Xe,V,8,T  Jq  Z 

where  12  is  a bounded  open  set  in  Rn  with  Lipschitz  boundary  (say  12  = (0,  256)2),  R 
is  the  rotation  matrix  in  terms  of  an  angle  9 , and  xe  is  the  characteristic  function  of 
E = {x  G 12|u(x)  > 0}.  This  minimization  is  over  all  the  characteristic  functions  of 
E in  BV(12)  (note  that  the  set  E varies  as  u evolves),  /t  G (0,  256],  9 G (— 7r,7r],  and 
T G [-256,256]  x [-256,256]. 

The  minimization  problem  (4.11)  involves  a weighted  total  variation  norm 
for  the  functions  with  finite  perimeters.  It  also  minimizes  more  than  one  argument. 
A special  case,  where  the  weight  function  is  one,  and  Xe  is  the  only  minimizing 
argument,  has  been  studied  in  [11].  To  study  the  existence  for  the  problem  (4.11),  it 
is  necessary  to  introduce  the  concept  of  weighted  total  variation  norms  for  functions 
of  bounded  variation.  Let  us  begin  with  recalling  the  definitions  of  functions  with 
bounded  variation  [29]. 

Definition  1:  Let  12  C Rn  be  an  open  set  and  let  f G L1^).  Define 


f | Vf | =:  sup{  f f(x)  div^(x)  dx}, 
Jn  <fre$  J n 
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where 

$ - {0  € Cj(D,Rn)|  |0(x)|  < 1,  on  ft}. 

Definition  2:  A function  f G L*(ft)  is  said  to  have  bounded  variation  in  ft,  if 
fQ  | Vf  | < oo.  We  define  BV(ft)  as  the  space  of  all  functions  in  L1  (^2)  with  bounded 
variation.  For  a function  f G BV(ft),  we  define 

||f||BV(n)  = ||f||Li(n)  + [ lVfl- 

Jn 

Definition  3:  A measurable  subset  E of  Rn  has  finite  perimeter  in  ft,  if  the  charac- 
teristic function  xe  G BV(ft). 

It  has  been  known  that  if  f G BV(ft),  then  for  any  t G R,  the  level  set 
Et  = {x  G ft|f(x)  > t}  has  finite  perimeter,  i.e.  XEt  G BV(ft). 

Next  we  define  the  weighted  total  variation  norm  with  the  weight  function  o(x). 
Definition  4:  Let  fl  C Rn  be  an  open  set.  Let  f G Lx(ft)  and  a(x)  be  positive 
valued  continuous  and  bounded  functions  on  ft.  The  weighted  total  variation  norm 
of  f with  the  weight  function  a(x),  denoted  by  a(x)|Vf|,  is  defined  by 

f ct(x)|Vf|  =:  sup  { f f(x)div</>(x)dx}, 

Jn  Jn 

where 

= {(j)  G Cj(ft,  Rn)|  |0(x)|  < o(x),  on  ft}.  (4-12) 

Lemma  5:  (see  [29])  Let  ft  C Rn  be  an  open  set  with  a Lipschitz  boundary.  Suppose 
that  fn  is  a bounded  sequence  in  BV(ft).  Then,  there  is  a subsequence  converging  in 
Lp(ft)  to  a function  f G BV(ft)  for  any  1 < p < 

Main  Theorem:  Let  ft  = (0,  256)2  and  let  C*  be  a differentiable  contour,  and  I G 
L°°(ft).  The  minimization  problem  (4.11)  has  a solution  xe  G BV(ft),  /i  G (0,256], 
6 G ( — 7r, 7r],  and  T G [-256,256]  x [-256,256]. 

Proof:  Let  Em  , /im,  Rm,  and  Tm  be  minimizing  sequences  of  (4.11).  Denote 


P(x)  = P(x,/x,R,T) 
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g(|VI|)(x)  + ^d2(MRx  + T). 


and 


Pm(x)  - P(x,  ^m,  Rm, Tm). 


Then,  as  m -4  oo 


(4.13) 


Since  g(|VI|)(x)  > l/(l  + K||I||L°°(n))  <,  for  a suitable  constant  K = K(/3,  a)  > 
0,  the  function  P is  bounded  below  by  a constants  p0  > 0.  Therefore,  y;Em  is  a bounded 
sequence  in  BV(Q).  By  Lemma  5,  there  exist  subsequences  of  XEm,  /hn,  #m>  and  Tm 
such  that,  without  changing  notation,  XEm  converges  to  a function  f in  L1^)  and 
a.e.  on  fb  Since  XEm  is  either  0 or  1,  f is  either  0 or  1 a.e..  We  may  view  f as  the 
characteristic  function  xe  of  a set  E.  Therefore,  we  have 


Moreover,  from  the  boundedness  of  //m,  Rm,  and  Tm,  there  exist  a constant 
/r  £ (0,256],  a constant  6 £ (— 7r,7r],  and  a constant  vector  T £ [0,256]  x [0,256], 
such  that 


Noticing  that  fx,  R,  and  T are  involved  in  the  continuous  function  P through  a distance 
function,  we  have  Pm  converges  to  P uniformly  on  fb  Therefore,  for  any  0 < e < p0 
there  exists  a constant  M = M(e)  > 0 such  that  if  m > M 


Then  by  the  definition  of  (ref.  (2))  we  have  $p_£  C <3>p.  Therefore,  for  any  fixed 
4>  £ $p_e,  if  m > M, 


XEm-tXE,  inL1^). 


/fin  ~ t ^ Tm  — * T. 


P - e < Pm  < P + e- 
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Figure  4.1:  I)  prior  shape  II)  image  I with  initial  contour  C0  III)  I with 

the  extracted  boundary  C.  IV)  C aligned  with  prior  shape  by  transformation 
M T). 


/ 


Pm|VXE 


Furthermore,  since  XEm  converges  in  Lx(0)  to  xe,  for  any  (f>  G $p- 

/ XEdiv^dx  = lim  / XEmdiv</> 

Jn  m^ooj  n 

The  combination  of  (4.14)  and  (4.15)  shows  for  any  (f>  £ <FP_£, 


Therefore, 


/ XEdiv</>dx  < lim  / Pm|VxEm|. 

Jn  Jn 


sup  / XEdiv0dx<  lim  / Pm|VxEm|- 
in 


By  letting  e — > 0,  and  then  using  the  definition  4 and  (4.11),  we  get 

/ P|DxeI  = SUP  / X p div./.dx 
Jn  Jn 


(4.14) 


(4.15) 


< lim 

m— KX) 


inf  P|Dxe|- 

XEeBV,^,R,T 


f Pm|VXEm 
Jn 

This  shows  that  solution  xe  £ BV(fi),  and  xe  Ab  and  T are  the  solutions  to  the 
minimization  problem  (4.11). 


67 


in 


Figure  4.2:  I)  High  Resolution  image  with  the  corpus  callosum  as  the  prior 
shape  II)  Two  low-resolution  images  R and  R with  the  initial  contour  (ellipse) 
and  the  extracted  boundaries  Ci  and  C2.  Ill)  Ci  and  C2  aligned  with  prior 
shape  S by  transformations  {pi,9i,Ti)  and  (//2,#2,T2). 

4.2.3  Numerical  Results 

We  have  tested  our  model  on  some  synthetic  images  and  functional  MR  brain 
images.  The  first  experiment(Fig.4.1)  we  present  here  is  to  detect  an  object  in  a 
binary  image  I containing  a half  disk.  The  image  is  given  in  (II)  and  the  shape  of 
the  object  to  be  detected  is  given  in  (I).  The  goal  is  to  find  a contour  C in  I that 
captures  high  gradient  (the  half  circle)  and  the  shape  prior  S.  Evolving  the  curve  Co, 
that  is  the  ellipse  in  (II),  according  to  the  system  of  equations  Eqs.(4.2)  - (4.5)  with 
the  parameters  A = 1,  po  = 1.4,#0  = 0,  T0  = (0,0),  we  get  the  boundary  C,  that 
is  the  contour  in  (III),  and  the  transformation  parameters  p = 1.25,  6 — 0.50,  and 
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I II 


Figure  4.3:  Functional  activation  maps  for  an  overt  word  generation  task 

obtained  by  (I)  AFNIs  intensity-based  registration  method,  (II)  the  proposed 
feature-based  method. 

T = (-2.14,  — 2.13)/128.  In  (IV),  the  transformed  edge  pRC  + T is  overlayed  on  the 
prior  shape  S.  From  (III),  we  can  see  that  the  contour  C captures  the  high  gradient 
in  the  image  I and  the  prior  shape  S,  even  though  the  gradient  information  is  totally 
missing  on  part  (the  line  segments)  of  the  contour  C. 

The  second  experiment  (Fig.4.2)  is  to  align  a set  of  time  series  fMR  brain 
images.  The  fMR  low-resolution  brain  images  were  acquired  on  a control  subject  in 
a 1.5  TGE  Sigma  scanner  with  a 2-spiral  gradient  echo  pulse  sequence  (TR/TE/FA 
=875ms/40ms/50deg,  FOV=180  mm,  matrix  128x128,  and  slice  thickness  of  5 mm). 
The  high-resolution  image  was  obtained  with  a 3D-SPGR  sequence  (TR/TE/FA 
=27ms  /7ms/45deg,  FOV=240mm,  matrix  256x256,  and  slice  thickness  of  1.3mm). 

The  shape  prior  S is  the  boundary  of  the  corpus  callosum,  shown  in  a high- 
resolution  brain  image  in  (I),  which  is  outlined  by  a radiologist.  We  applied  the 
proposed  model  to  find  the  boundary  of  the  poorly  defined  corpus  callosum  in  each 
of  the  time  series  images  from  the  same  subject.  In  addition  to  this,  the  spatial  trans- 
form that  maps  the  new  segmented  contour  to  the  given  contour  S is  also  obtained. 
This  information  leads  to  the  co-registration  of  the  time  series  images.  In  (II),  two 
temporally  adjacent  time  series  fMR  brain  images  from  the  same  subject  are  shown. 
The  ellipses  in  these  two  images  are  the  initial  contours.  Evolving  them  according  to 
Eqs.(4.2)  - (4.5)  with  initial  values,  po  = 1-5,  #o  = 0,  To  = (0,0),  we  obtained  the 
boundaries  in  (II) (the  curves  within  the  ellipses),  denoted  by  Ci  and  C2,  respectively. 
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Parts  of  them  are  captured  by  the  intensity  gradient.  The  corresponding  transforma- 
tion parameters  are  /xi  = 1.3585, 0i  = 0.1169,  Ti  = (37.0268, 44. 8614) /128  that  maps 
Ci  to  S,  and  /x2  = 1.38  1 5, 02  = 0.1027,  T2  = (36.1594, 42.8034) /128  that  maps  C2  to  S. 
The  contours  in  (III)  are  /xiRiCi  -t-Ti,  and  /i2R2C2 + T2,  respectively.  By  comparing 
the  figures  in  (III)  with  (I),  we  can  see  that  the  contours  Q (i  =1,2)  are  captured  by 
the  shape  prior  S.  Finally,  the  transformation  C2  = 1Ci+/i2  *R2  X(Ti—  T2) 

provides  the  alignment  of  these  two  time  series  images. 

Fig(4.3)  shows  the  functional  activation  maps  for  an  overt  word  generation 
task  obtained  by  (I)  using  AFNIs  intensity-based  registration  method  and  (II)  using 
our  feature-based  method.  In  this  experiment,  the  task  was  single-word  generation 
to  a series  of  randomly  spaced  category  cues.  Four  functional  runs  were  acquired 
with  150  images  in  each  run,  and  1.5  sec. per  image.  The  red  dots  in  the  figures  rep- 
resent high  signal.  In  Fig(4.3),  (I)  and  (II)  show  that  the  feature-based  registration 
method  produces  a significant  decrease  in  the  motion-related  noise  compared  to  the 
AFNI  technique.  Regional  activation  in  medial  cortex  at  and  near  the  cross  hair 
results  from  local  blood  flow  and  oxygenation  changes  which  affect  the  MR  signal 
intensity  subsequent  to  neuronal  activity.  These  correspond  to  areas  also  observed  in 
covert  language  generation  studies.  Compared  to  AFNI,  the  feature-based  registra- 
tion method  reduces  motion  artifacts,  which  are  predominant  in  areas  of  high  signal 
contrast. 


The  M-S  functional  is  one  of  most  widely  studied  mathematical  models  in  im- 
age processing  and  addresses  the  problems  of  image  reconstruction  and  segmentation, 
simultaneously.  It  is  formulated  as  follows: 


4.3.  M-S  functional  Based  Segmentation  with  Shape 
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where, 

C denotes  the  closed  segmented  curve,  I denotes  the  given  image  and  f is  the  piecewise 
approximation  to  I,  with  discontinuities  only  along  C. 

Let  R and  Rc  denote  the  regions  inside  and  outside  C,  respectively.  We  can 
find  the  required  piecewise  smooth  f by  computing  two  smooth  functions  fo  and  fc 
such  that, 


{fo  in  R 

(4.16) 

fc  in  Rc 

Now,  the  [M-S]  functional  can  be  re-written  to  simplify  computations: 

E(f0,fc,C)  = 

[ (f0  — I)2dx  + [ (fc  — I)2dx  4-  a [ |Vf0|2dx  + a [ |Vfc|2d x + 7/  ds  (4.17) 
Jr  J rc  Jr  J rc  J c 

Like  most  region  based  approaches,  the  [M-S]  segmentation  model  enjoys 
greater  robustness  to  noise(by  avoiding  derivatives  of  I)  and  initial  contour  place- 
ment (use  of  global  image-intensity  information)  than  edge-based  methods  like  the 
GAC. 

But  the  model  operates  with  an  implicit  assumption  that  the  boundaries  of  I separate 
homogeneous  regions,  which  is  not  true  in  the  case  of  fMR  images.  We  can  add  the 
shape  term  to  energy  (4.17)  to  alleviate  this  issue: 

E(f0,fc,C,/qR,T)  = [ (f0  - I)2dx  + [ (fc  — I)2dx  + a [ |Vf0|2dx 
Jr  J rc  Jr 

+a  ( |Vfc|2dx  + j (f  d2(pRC  + T)  ds  (4.18) 

J Rc  J C 

In  order  to  allow  C to  change  its  topology  and  also  to  re-write  the  above 
functional  as  a single  integral  over  fi,  we  use  a level-set  formulation  as  in  ([11]). 

A contour  C is  represented  by  the  zero  level  set  of  a lipschitz  function  u such 
that  Rc  = (x|u(x)  > 0},  the  set  outside  C.  Now,  we  can  write  the  energy  in  Eq.(4.18) 


as 
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E(f0,fc,u,/x,R,T)=  [ (1  — H(u))(f0  — I)2dx  + [ H(u)(fc  - I)2dx 
Jn  Jn 

+a  f (1  — H(u))|Vfo|2dx  + a [ H(u)|Vfc|2dx  + 7 [ <5(u)d2(//RC  + T)|Vu|  d)(4.19) 
*/  r 2 «/  n n 

where,  H(z)  is  the  heaviside  function  and  6 = H'  in  the  distributional  sense. 

4.3.1  Minimization  of  Energy 

The  above  minimization  problem  is  solved  by  looking  at  the  steady  state  of 
the  following  system: 

= -<5u[F  - 7div(d2|^j)],  x G Q,  t > 0 (4.20) 

— ^ = 0,  x G dQ,  t > 0, 
an 

u(x,  0)  = u0(x),x  G and 
F = [fa  - I)2  - (fc  - I)2]  + a[\ Vf0|2  + |Vfc|2]. 


or 

= -[(1  - H(u))(f0  - I)  - adiv((l  - H(u))Vf0)),  x6fl,t>0  (4.21) 

or 

-7-^  = 0,  x E C,  t > 0, 
an 

f0(x,  0)  = f0(x),x  G R. 


= -[H(u)(fc  - I)  - adiv(H(u)Vfc)),  xEfi,  t>0 

Of 

— ^ - = 0,  x G C,  t > 0, 
an 

fc(x,  0)  = fc(x),x  G Rc. 


(4.22) 
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/it  = — A / d(u)|Vu|dVd  • Rx  dx,  /i(0)  = /xq, 


(4.23) 


n 


(4.24) 


Tt  = -A  [ <5(u)|Vu|dVd  dx,  T(0)  = T0, 


(4.25) 


Jn 


where  n is  the  outward  unit  normal  to  C. 

To  numerically  solve  the  above  equations,  we  use  an  explicit  finite  difference 
scheme.  In  our  implementation,  we  replace  H and  d in  Eq(4.20)-Eq(4.25)  by  the 
slightly  regularized  versions  of  them,  denoted  as  Hf  and  5e  respectively(as  in  [11]). 

4.3.2  Numerical  Results 

Similar  experiments  as  for  the  GAC  model  were  conducted  on  synthetic  and 
fMR  brain  images.  The  first  experiment(Fig.4.4)  we  present  here  is  to  detect  an 
object  in  a binary  image  I(x,y)  containing  a disk,  partially  occluded  by  a square  (I). 
Our  goal  is  to  find  a contour  C in  (I)  that  captures  the  edge  information  and  is  similar 
to  the  shape  prior  S(II). Evolving  an  initial  curve  C0(ellipse  in  (I))  according  to  the 
system  of  equations  Eq(4.20)-Eq(4.25)  we  get  the  solution  C(curve  in  (III) )that  has 
completed  the  missing  portion  of  the  disc  using  the  shape  information. 

The  second  experiment,  Fig.  4.5,  is  to  align  a set  of  time  series  fMR  brain 
images.  We  use  two  temporally  adjacent  time  series  fMR  brain  images  from  the 
same  activation  experiment  that  was  outlined  in  the  previous  section  on  the  GAC 
model.  The  ellipses  in  (I)  are  the  initial  contours.  Evolving  them  according  to 
Eq(4.20)-Eq(4.25),  we  find  the  contours  in  (II),  for  each  image,  denoted  by  Cl  and 
C2,  respectively.  The  alignment  of  the  transformed  contours  /r, R;Cj  + T;  with  the 
shape  S is  shown  in  (III). 
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B 

I II  ill 

Figure  4.4:  Filling  missing  boundary  with  prior  shape  (I)  Image  I with  initial 
contour  (II)  Prior  shape  (III)  Extracted  boundary  C. 

4.4.  Conclusion 

We  proposed  the  addition  of  a shape  prior  to  an  active  contour  in  a variational 
framework.  The  key  idea  was  to  introduce  an  energy  term  which  measured  the 
closeness  of  the  active  contour  to  a prior  contour  in  a way  that  is  independent  of 
the  location,  orientation,  and  scale  of  the  contour.  In  experiments  , the  shape-based 
active  contour  could  reliably  segment  images  in  which  the  complete  boundary  was 
either  missing  or  was  low-resolution  and  low-contrast.  Besides  the  segmentation, 
the  algorithm  also  provides  estimates  of  translation,  rotation,  scale  that  map  the 
active  contour  to  the  prior  shape.  These  estimates  are  useful  in  aligning  images.  The 
extension  to  3-d  images  can  be  done  (the  shapes  are  aligned  using  3-d  translation, 
rotation,  and  scaling)  and  will  be  the  subject  of  future  reports. 
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Prior  shape  S 


Figure  4.5:  Registering  2 fMR  images  (I)  Image  with  initial  contour  (II) 

Extracted  boundaries  Ci  and  C2  (III)  C\  aligned  with  shape  S in  the  High- 
Resolution  image. 
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